DTIC  COPY 


AFRL-RV-H  A-TR-2007-1 1 52 


Equator  and  High-Latitude  Ionosphere-to-Magnetosphere 
Research 


B.  W.  Reinisch 
G.  S.  Sales 
V.  Paznukhov 
I.  A.  Galkin 
D.  F.  Altadill 
G.  Khmyrov 


Center  for  Atmospheric  Research 
University  of  Massachusetts  Lowell 
600  Suffolk  Street 
Lowell,  MA  01854 


Scientific  Report  No.  1 
30  October  2007 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED. 


AIR  FORCE  RESEARCH  LABORATORY 


Space  Vehicles  Directorate 
29  Randolph  Road 

AIR  FORCE  MATERIEL  COMMAND 
Hanscom  AFB,  MA  01731-3010 


20080624  204 


NOTICE  AND  SIGNATURE  PAGE 


Using  Government  drawings,  specifications,  or  other  data  included  in  this  document  for 
any  purpose  other  than  Government  procurement  does  not  in  any  way  obligate  the  U  S. 
Government.  The  fact  that  the  Government  formulated  or  supplied  the  drawings, 
specifications,  or  other  data  does  not  license  the  holder  or  any  other  person  or  corporation; 
or  convey  any  rights  or  permission  to  manufacture,  use,  or  sell  any  patented  invention  that 
may  relate  to  them. 

This  report  was  cleared  for  public  release  and  is  available  to  the  general  public,  including 
foreign  nationals.  Qualified  requestors  may  obtain  additional  copies  from  the  Defense 
Technical  Information  Center  (DTIC)  (http://www.dtic.min.  All  others  should  apply  to  the 
National  Technical  Information  Service. 


AFRL-RV-HA-TR-2007-1 152  HAS  BEEN  REVIEWED  AND  IS  APPROVED  FOR 
PUBLICATION  IN  ACCORDANCE  WITH  ASSIGNED  DISTRIBUTION  STATEMENT. 


/  /  Signature/  /  /  /Signature/  / 


KENNETH  R.  WALKER  JOEL  MOZER,  Chief 

Contract  Manager  Space  Weather  Center  of  Excellence 


This  report  is  published  in  the  interest  of  scientific  and  technical  information  exchange,  and 
its  publication  does  not  constitute  the  Government’s  approval  or  disapproval  of  its  ideas  or 
findings. 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

Hcreportjngburdenfor  the  collection  of  information  ,s  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  search mg  existing  data  sources,  gathenng  and  maintaining 

i  needed^  HHS.  rev'e™n9thei  coJectK>n  of  information  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information  including  suggestions  for  red u< 

burden  to  Department  of  Defense.  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188).  1215  Jefferson  Davis  Highway  Suite  1204  Adu^toTvA  222(7 

PLEASE  TO  NOTRFTURN'YOUR^FORWrro  THE* ABOVE  ADDRESS.*  *  SU^to  «  «“ «"■  <-p,y  wm,  a  cCec^oo  a.  iirfoonaton  if  it  does  not  a  cuo 

REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

1-10-2007  Scientific  Report  No.  1 

3.  DATES  COVERED  (From  -  To) 

25  Aug  2006  -  25  Aug  2007 

riTLE  AND  SUBTITLE  - 

[uator  and  High- Latitude  Ionosphere-to-Magnetosphere  Research 

5a.  CONTRACT  NUMBER 

FA8718-06-C-0072 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

62601F 

AUTHOR(S)  ~  - - 

W.  Reinisch,  G.S.  Sales,  V.  Paznukhov,  I. A.  Galkin,  D.F.  Altadill,  G.  Khmyrov 

5d.  PROJECT  NUMBER 

4827 

5e.  TASK  NUMBER 

HR 

5f.  WORK  UNIT  NUMBER 

A1 

PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES)  ~ 

inter  for  Atmospheric  Research 
tiversity  of  Massachusetts  Lowell 

10  Suffolk  Street 
>well,  MA  01854 

8.  PERFORMING  ORGANIZATION  REPOf 

NUMBER 

SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESSES)  " 

.r  Force  Research  Laboratory 

)  Randolph  Road 

tnscom  AFB,  MA  01731-3010 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

AFRL/RVBXI 

11.  SPONSOR/MONITOR’S  REPORT 

NUMBER(S) 

AFRL-RV-HA-TR-2 007- 1152 

>proved  for  public  release;  distribution  unlimited. 


SUPPLEMENTARY  NOTES 


.  ABSTRACT  - 

is  report  describes  significant  improvements  in  digital  sounder  technology  as  applied  to  ionospher: 
search  in  the  equatorial,  mid  and  high-latitude  regions.  These  new  developments  are  in  the  areas  < 
*no spheric  drift,  radio-wave  absorption,  precision  group-height  and  the  specification  of  sounder 
asurement  uncertainty.  Ultimately  these  techniques  are  to  be  included  into  routine  sounder 
►erations  and  ARTIST  performance.  These  developments  required  research,  which  is  reported  here,  int- 
ie  structure  of  the  ionospheric  layers,  including  the  D- region,  and  the  E  and  F  layers.  In  addition 
jor  effort  under  this  project  involved  ionospheric  profile  validation;  that  is  comparing  sounder 
asured  electron  density  profiles  with  the  results  of  other  techniques.  These  comparisons  were  mad< 
ing  the  data  from  other  researchers  including  tomographic,  occultation  and  DMSP  UV  measurements, 
nally,  the  effects  of  strong  magnetic  storms  on  a  dense  network  of  sounders  in  Europe  are  presenter 
rge  changes  in  the  ionosphere  were  observed  in  the  daytime  on  disturbed  days  and  the  latitudinal 
sponse  of  ionosphere  was  analyzed  and  compared  to  similar  observations  at  American  longitudes  wher 
ry  different  response  was  observed. 

SUBJECT  TERMS  - - - 

ignetic  storms.  Ionospheric  drift.  Radio  absorption.  Profile  validation,  ARTIST  performam 


SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT 

18.  NUMBER 

OF  PAGES 

19a.  NAME  OF  RESPONSIBLE  PERS 

Kenneth  R.  Walker 

REPORT 

U 

b.  ABSTRACT 

U 

c.  THIS  PAGE 

U 

SAR 

76 

19b.  TELEPHONE  NUMBER  (Include  t 
code) 

Standard  Form  298  (Rev.  8-38) 

Prescribed  by  ANSI  Std.  239.18 


Contents 


1.  INTRODUCTION  1 

2.  COMPARISON  OF  THE  IONOSPHERIC  PROFILES  DERIVED  USING  A  RADIO  1 
OCCULT ATION  TECHNIQUE  TO  THE  DIGISONDE  DATA  1 

2.1  Introduction  1 

2.2  Data  Collection  2 

2.3  Data  Analysis  3 

2.4  Conclusions  and  Future  Efforts  9 

3.  ARTIST  UNCERTAINTY  SPECIFICATION  1 0 

3.1  Introduction  10 

3.2  Approach  10 

3.3  Statistical  Analysis  of  ARTIST  Performance  Against  Manual  Data  1 1 

3.3.1  General  Approach  11 

3.3.3  ARTIST  Confidence  Level  of  Autoscaling  and  Uncertainty  12 

3.3.4  Calculating  Uncertainty  Bounds  13 

3.3.5  Boulder  Digisonde  Dataset  14 

3.3.6  Classification  by  Spread  F  Severity  Level  and  CLA  1 5 

3.3.7  Uncertainty  Associated  with  Different  Versions  of  ARTIST  1 7 

3.3.8  Station-Specific  Uncertainty  18 

3.4  Reporting  Uncertainties  in  SAO  Format  19 

3.4.1  Uncertainty  Bounds  for  N(h)  Profiles  19 

3.4.2  Software  Modifications  19 

3.5  Future  Efforts  20 

4.  DIGISONDE  GROUND-TRUTH  DATA  FOR  THE  NWRA  TOMOGRAPHY  2 1 

VALIDATION  PROJECT 

5.  PRECISION  GROUP  HEIGHT  MEASUREMENTS  OF  E  LAYER  VIRTUAL  2 1 

HEIGHTS 

5.1  Introduction  21 

5.2  The  Virtual  Height  Phase  Technique  22 

5.3  Measuring  E  :Layer  Heights  at  Millstone  Hill  25 

6.  PROCESSING  OF  DIGISONDE  DRIFT  MEASUREMENTS  29 

6.1  Introduction  29 

6.2  Data  Analysis  30 

6.3  Future  Efforts  on  Drift  Data  Analysis  34 

6.4  DDA  Software  Upgrade  35 

7.  STUDY  OF  THE  IONOSPHERIC  RESPONSE  TO  STRONG  GEOMAGNETIC  35 
STORMS 

7.1  Introduction  35 

7.2  Data  Analysis  36 

iii 


7.2  Data  analysis  36 

7.3  Discussion  42 

7.4  Summary  44 

8.  ROUTINE  HF  ABSORPTION  MEASUREMENTS  USING  DIGISONDES  45 

8.1  Approach  45 

8.2  System  calibration  46 

9.  FOCUSING  AND  REFLECTION  CALCULATIONS  48 

9. 1  Focusing  gain  48 

9.2  Nighttime  F-layer  reflectivity  52 

9.3  Reflectivity  theory  54 

1 0.  CALIBRATION  /  VALIDATION  OF  DMSP  UV  MEASUREMENTS  OF  THE  F2  57 

PEAK  CHARACTERISTICS  WITH  DIGISONDE  MEASUREMENTS 

11.  PUBLICATIONS  61 

REFERENCES  63 


IV 


Figures 


1 .  A  COSMIC  satellite  pass  over  the  Athens  digisonde  on  December  21,  2006.  The  blue  3 
symbols  show  radial  projections  of  the  ray  path  tangent  points.  The  red  star  indicates  the  station 
location. 

2.  Athens  ionograms  recorded  during  the  occultation  on  December  2 1 ,  2006  at  0800  UT  4 

and  0815  UT. 

3.  Comparison  of  the  digisonde  and  COSMIC  RO  profiles  for  Athens.  The  plots  show  the  5 
electron  density  (top  figure)  and  plasma  frequency  (bottom)  as  a  function  of  height. 

4.  Comparison  of  the  Ascension  Island  profiles  derived  from  digisonde  and  COSMIC  RO  6 
measurements  (plots  on  the  left)  obtained  on  December  20, 2006.  On  the  right  is  the 
measured  ionogram  demonstrating  spread  F  conditions. 

5.  Comparison  of  the  Kwajalein  Island  profiles  obtained  from  digisonde  and  COSMIC  RO  7 
measurements  on  December  21,  2006. 

6.  RO  tangent  points  (blue)  near  Zhigansk  and  Yakutsk  digisonde  stations  on  8 

December  2 1 , 2006 

7.  RO  measurements  on  December  2 1 ,  2006  in  the  vicinity  of  Zhigansk  and  Y akutsk  8 

digisonde  stations. 


8.  Comparison  of  foF2  values  deduced  from  simultaneous  digisonde  and  RO  measurements.  9 

9.  Point-by-point  uncertainty  of  the  ARTIST  electron  density  profile  is  obtained  by  10 

Calculating  inner  and  outer  boundaries  enclosing  the  profile.  Both  boundaries  rest  on  five  anchor 
points  whose  uncertainties  are  known  from  statistical  analysis  of  manually  evaluated  differences 
between  automatically  and  manually  scaled  values. 

10.  A  digisonde  ionogram  recorded  during  severe  spread  F  conditions  in  Jicamarca,  Peru  (left).  12 
If  multiple  echoes  are  resolved  and  presented  as  individual  edgels  (right),  the  trace  extraction 
algorithm  is  overwhelmed  with  numerous  possibilities  of  grouping  edgels  to  traces. 

1 1 .  Cumulative  difference  function  of  the  absolute  |autoscaled-manual|  differences  for  a  rapid  1 3 
evaluation  of  ARTIST  uncertainty.  ARTIST-5  reports  of  foF2  for  Grahamstown  have  0. 1 5  MHz 
uncertainty  at  95%  level. 
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12.  Uncertainty  bounds  for  foF2  derived  from  the  frequency  histogram  of  (ARTIST-manual)  14 
differences  for  Boulder  Digisonde  256  ionograms.  For  12,612  ionograms  scaled  with  confidence 
(87%  of  total  14,501  ionograms),  the  uncertainty  bounds  at  5%  percentile  are 

-0.25  and  +0.35  MHz. 

13.  Availability  of  manually  interpreted  data  for  the  Boulder  Digisonde  station.  The  plot  1 5 
shows  foF2  for  every  processed  ionogram  for  years  2004-2006. 

14.  ARTIST  4.5  uncertainty  calculated  for  the  winter  and  summer  time  for  Boulder  digisonde.  16 
Plots  show  foF2  cumulative  differences  calculated  for  nighttime  (black)  and  daytime  (red) 
ionograms  for  each  period. 

1 5.  Uncertainty  of  foF2  autoscaled  value  for  Gakona,  AL  digisonde  as  a  function  of  1 6 

automatically  detected  level  of  spread  F  conditions.  Comparison  of  ARTIST  5.0.2-b7  results  for 
5474  manually  scaled  ionograms  during  the  2002-2006  period. 

16.  Cumulative  difference  characteristics  for  Boulder  2004  ionograms.  This  figure  compares  1 7 
the  results  obtained  with  ARTIST5.0  left)  and  ARTIST4.5  (right).  Numerical  values  for  foF2, 
foFl,  and  foE  scaling  uncertainties  (calculated  at  95%  level)  are  listed  as  well. 

1 7.  Calculated  profile  uncertainty  boundaries  shown  in  the  SAO  Explorer.  20 

1 8.  Distribution  of  virtual  heights  during  one  sounding  made  on  14  September  2006  at  26 
1422  UT.  A  majority  of  the  calculated  heights  are  concentrated  within  a  single  statistical  range 
bin  at  107.5±0.5  km. 

19.  Daily  E  region  precision  group  heights  at  Millstone  Hill  as  function  of  time.  The  27 

horizontal  axes  represent  time  (Universal  Time),  and  the  vertical  axes  represent  virtual 

height  (km).  Dashed  circles  highlight  significant  decreases  in  E  layer  virtual  heights  and 
the  black  arrows  point  to  hook-shaped  disturbances. 

20.  Average  daily  variations  of  the  h’E  obtained  from  the  20  days  of  September.  Vertical  28 
error  bars  indicate  the  standard  deviation  obtained  at  a  given  time. 

21 .  Daily  pattern  of  the  vertical  velocity  data  from  the  Ebro  digisonde  for  December  2004.  30 

The  15  minute  median  values  of  Vz  are  shown  as  gray  dots.  The  vertical  error  bars  indicate  the 
range  of  80%  of  Vz  values  for  each  given  time.  The  solid  line  represents  the  daily  pattern 
calculated  as  a  sum  of  the  primary  diurnal  harmonics.  The  vertical  dashed  lines  indicate  the 
sunset  (SS)  and  sunrise  (SR)  above  the  station  at  an  altitude  of  300  km. 

22.  Seasonal  variations  of  the  amplitudes  (left)  and  phases  (right)  of  the  daily  pattern  of  the  3 1 
Vz  in  the  F  region.  The  top  plots  show  the  diurnal  component,  the  middle  and  bottom  ones 
correspond  to  the  semidiurnal  and  terdiumal  component,  respectively. 

23.  The  left  panel  shows  the  seasonal  variation  of  the  phase  of  the  semidiurnal  harmonic  3 1 
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(full  line  with  full  dots),  the  time  of  the  minimum  values  of  the  daily  pattern  of  Vz  (dashed  line 
with  open  squares)  and  the  sunrise  time  at  300  km  height  (full  line  with  open  circles).  The  right 
panel  shows  the  scatter  plot  of  the  semidiurnal  phase  against  the  sunrise  time  for  a  given  month. 
Solid  line  is  the  best  linear  fit  with  the  correlation  coefficient  indicated. 

24.  Monthly  averaged  east-west  drift  velocity  data  measured  at  Jicamarca  station  in  2006.  33 

25.  Amplitude  spectra  of  the  east-west  plasma  velocity  measured  at  Jicamarca  in  2006.  34 

Bottom  axis  shows  the  frequency  in  pHz  while  the  top  axis  shows  the  corresponding  time 
periods  in  hours. 

26.  Geomagnetic  storm  of  24  August  2005.  The  top  panel  shows  SYM-H  index  variations,  37 
and  the  next  panels  show  the  solar  wind  plasma  density,  plasma  speed,  and  the  last  two 

panels  show  thw  interplanetary  magnetic  field  components.  The  solar  wind  and  IMF  data 
were  taken  using  the  ACE  satellite  and  were  shifted  by  40  minutes  to  account  for  the  position  of 
the  satellite. 

27.  Electron  density  profile  variations  as  a  function  of  time  for  the  listed  European  stations  38 
(see  Table  5).  The  right  panel  shows  storm  day  (August  24,  2005)  data,  while  the  left  panel 
shows  the  quiet  day  average  data.  Electron  density  is  presented  in  terms  of  equivalent  plasma 
frequency  (color  bar).  Sunrise  and  sunset  times  are  indicated  for  the  Ebro  station  with  red  and 
black  circles  correspondingly. 

28.  Variation  of  foF2  and  hmF2  layer  parameters  on  24  August  2005  compared  to  the  quiet  39 
day  averages  (black  curves).  Top  panel  shows  DST  index. 

29.  American  sector  data  for  24  August  2005  storm.  Plots  show  variations  of  foF2  (left)  and  40 
hmF2  (right)  as  compared  to  the  quiet  day  averages  (black  curves). 

30.  Difference  between  hmF2  values  measured  on  the  event  day  and  the  quite  day  average  as  41 
a  function  of  time  after  storm  commencement. 

3 1 .  August  24,  2005  storm.  Color  scale  shows  the  energy  input  determined  from  radar  and  43 
DMSP  satellite  measurement.  These  data  are  taken  from  OVATION  database.  The 

locations  of  selected  digisonde  stations  as  well  as  Bjomoya  magnetic  observatory  are 
shown  in  the  plot. 

32.  Variations  of  the  AhmF2  (black)  and  AfoF2  (blue)  parameters  measured  in  the  European  44 
sector.  The  increase  in  the  electron  density  lags  height  increases  at  all  stations.  A  vertical  line 
indicates  the  storm  commencement  time  (SST).  The  sloping  lines  connect  the  starts  of  the 
increases  in  hmF2  (black  line)  and  foF2  (blue  line). 

33.  Number  of  amplitude  measurements  used  to  determine  the  average  PTGTGR  as  a  function  46 
of  sounding  frequency. 
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34.  Histogram  of  the  digisonde  derived  PTGTGR  values  over  the  nighttime  period  from  May  47 
through  August  1005  for  2.9  MHz. 

35.  Illustration  of  focusing  and  defocusing  gain  associated  with  an  irregular  reflecting  surface.  49 

36.  Distribution  of  100  realizations  of  an  irregular  reflecting  surface  with  RMS  amplitudes  50 
varying  from  0.1  through  5.0  km.  In  all  cases  the  horizontal  scale  size  was  »  6.8  km. 

37.  Percentage  positive  and  negative  focusing  gain  as  a  function  of  horizontal  structure  size  51 
for  several  vertical  RMS  amplitudes:  a)  5  km,  b)  10  km,  c)  15  km  and  d)  20  km. 

38.  Boulder  measured  nighttime  loss  (05-09  UT)  on  five  sequential  days.  53 

39.  Mean  loss  over  the  hours  from  05-09  UT  on  sequential  nights  from  3  May  -  7  May  2005  54 
(dots).  Calculated  mean  F-layer  slope  (km/MHz)  using  the  ionograms  for  the  period  from  1 
May-7  May  2005  (squares). 

40.  Reflection  index  ( 20  log  \R\ )  as  a  function  of  sounding  frequency  for  normalized  layer  55 

thicknesses  varying  from  2  (thinnest)  to  250  (thickest)  in  three  factors  of  5.  (Rawer  and  Suchy,  p. 
194,  Figure  64) 

41 .  Reflection  loss  (dB)  as  a  function  of  normalized  layer  thickness  for  two  collision  56 

frequencies. 
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1.  INTRODUCTION 


This  research  contributes  in  critical  areas  toward  the  goals  of  (1)  ionospheric  effects  on  DoD 
systems  research  and  (2)  ionospheric  research  technology. 

Based  on  the  proposal  submitted  to  the  Air  Force,  this  year’s  efforts,  as  described  in  this 
report,  involve  support  to  the  ionospheric  specification  objectives  of  the  Air  Force  Research 
Laboratory  (AFRL).  The  University  of  Massachusetts  Lowell  Center  for  Atmospheric  Research 
(UMLCAR)  has  taken  an  approach  that  addresses  the  specification  of  ionospheric  parameters  on 
a  global  scale;  a  goal  that  is  particularly  facilitated  using  the  digisonde  system  pioneered  by 
UMLCAR.  Global  ionospheric  modeling  is  a  major  part  of  space  weather  forecasting  and  global 
communications  progress,  and  our  support  for  these  goals  is  presented  here.  Use  of  the 
ubiquitous  digisonde  offers  the  best  tool  for  real  time  ionospheric  assimilative  modeling  as  well 
as  support  to  other  systems  that  require  verification  and  validation.  A  significant  part  of  our 
research  effort  was  providing  the  necessary  validation  of  other  methods  of  ionospheric 
specification  by  using  the  digisonde  measurements  as  the  “truth”  against  which  the  performance 
of  other  systems  was  compared.  These  include  cooperation  with  groups  making  ionospheric 
radio  occultation  and  tomography,  and  UV  measurements  (Sections  2, 4,  and  10).  The  Center  is 
also  evaluating  the  digisonde  drift  measurements  against  incoherent  scatter  radar  plasma  drift 
measurements  (Section  6).  The  second  major  research  concentration  in  this  report  involves  the 
development  of  techniques  that  expand  global  ionospheric  specification  (Sections  3,  5,  and  8). 
These  cover  the  areas  of  digisonde  profile  uncertainty,  improved  ionogram  virtual  height 
measurements,  and  using  digisondes  to  carry  out  routine  measurements  of  the  absorption  of  HF 
radio  waves  with  the  aim  of  predicting  system  outages  with  improved  sensitvity. 

Finally,  this  report  presents  the  study  of  the  effects  of  geomagnetic  storms  on  the  structure  of 
the  ionosphere  over  a  wide  range  of  latitudes  and  longitudes  (Section  7),  again,  addressing  the 
goals  of  global  modeling. 


2.  COMPARISON  OF  THE  IONOSPHERIC  PROFILES  DERIVED  USING  A  RADIO 
OCCULTATION  TECHNIQUE  TO  THE  DIGISONDE  DATA 

2.1  Introduction 

The  development  of  global  ionospheric  models  in  the  frame  of  the  “space  weather”  concept 
presents  the  challenge  of  establishing  a  global  observing  network  for  monitoring  the  Earth’s 
ionosphere.  One  promising  technique  suitable  for  establishing  a  foundation  for  such  an  ambitious 
project  is  the  radio  occultation  (RO)  method  that  allows  reconstructing  ionospheric  density 
profiles  over  a  large  altitude  range  and  has  the  natural  advantage  of  making  measurements  on  a 
global  scale.  Of  course,  a  significant  number  of  satellites  is  needed  to  provide  adequate  coverage 
of  the  global  ionosphere  but,  with  advances  in  space  technology,  this  could  eventually  be 
achievable  at  a  reasonable  cost.  At  the  present  stage  of  development  and  testing  of  the  RO 
method  for  the  ionosphere  sensing,  it  is  most  important  to  verify  the  accuracy  of  these 
measurements  using  an  established  ionospheric  technique.  Over  the  past  year,  we  have  worked 
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on  the  verification  of  the  RO  measurements  made  on  the  COSMIC  satellites  by  using  ground- 
based  digisonde  profile  measurements  at  various  locations.  Our  results  are  reported  here. 

2.2  Data  Collection 

The  COSMIC/FORMOSAT-3  (Constellation  Observing  System  for  Meteorology, 

Ionosphere,  and  Climate,  and  Taiwan's  Formosa  Satellite  Mission  #3)  is  a  joint  Taiwan-U.S 
program.  The  project  was  launched  in  December  of  2005,  and  is  expected  to  continue  providing 
data  for  five  years.  The  COSMIC/FORMOSAT-3  program  consists  of  the  six  spacecraft,  each 
with  three  instruments,  including  a  GPS  RO  receiver,  an  ionospheric  photometer,  and  a  tri-band 
beacon.  The  mission  is  aimed  at  space  weather  and  climate  research  and  forecasting,  as  well  as 
geodesy  and  gravity  research. 

On  20-21  December  2006,  a  joint  multi-instrument  campaign  was  carried  out  in  support  of 
the  COSMIC/FORMOSAT  mission.  One  of  the  primary  objectives  of  this  campaign  was  to 
validate  the  electron  density  profiles  determined  using  the  radio  occultation  technique.  Among 
the  instruments  participating  in  the  campaign  were  40  digisondes  distributed  around  the  globe. 
From  December  19  until  December  22,  2006,  most  of  the  digisonde  stations  in  this  program 
increased  their  measurement  cadence  to  5  minutes.  All  collected  data  were  archived  in  the 
Digital  Ionogram  Database  (http://ulcar.uml.edu/DIDBase/),  and  the  World  Data  Centers 
archived  the  auto-scaled  characteristics  and  profiles,  making  these  data  available  to  the  scientific 
community. 

Our  aim  was  to  establish  under  what  conditions  the  RO  technique  is  able  to  correctly  derive 
the  F2  layer  profiles,  especially,  the  peak  characteristics  foF2  and  hmF2,  as  well  as  the  E  layer 
characteristics  foE  and  hmE,  and  what  are  the  typical  errors.  In  cases  of  good  agreement  of  the 
F2  peak  characteristics,  the  shapes  of  the  bottom  and  topside  profiles  were  compared.  We  started 
the  comparison  by  preparing  a  list  of  times  where  RO  profiles  over  or  near  Digisonde  stations 
were  expected.  For  these  times,  the  digisonde  autoscalings  were  manually  edited  and  the  edited 
data  were  added  to  the  DIDBase  archive.  We  have  processed  about  70%  of  the  predicted  cases 
for  the  December  20-21,  2006,  campaign.  Unfortunately,  not  all  the  predicted  satellite 
occultations  produced  successful  electron  density  profiles,  especially  at  high  latitude  and 
equatorial  regions.  Thus,  we  were  able  to  collect  only  25  simultaneous  measurements  for  which 
both  digisonde  and  RO  profiles  existed. 

To  be  exact,  the  RO  density  profile  does  not  represent  an  actual  vertical  profile  for  the 
tangent  point,  but  rather  an  average  density  profile  representative  for  the  ray  path  tangent  points. 
The  assumption  of  local  spherical  symmetry  of  the  density  distribution  in  a  large  region  (up  to  a 
few  thousand  kilometers)  is  used  in  order  to  retrieve  the  vertical  profiles.  The  size  of  the  region 
for  which  this  assumption  was  applied  is  characterized  by  a  smear  parameter,  which  is  the 
horizontal  distance  between  the  top  and  bottom  tangent  points  of  the  measurements.  In  the 
December  2006  campaign,  for  the  observations  of  interest,  smear  factors  varied  from  several 
hundred  to  a  few  thousand  kilometers.  The  next  section  summarizes  the  results  of  the 
comparisons. 
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2.3  Data  Analysis 


The  first  example  used  the  occultation  measurements  from  0758  UT  to  0815  UT  on  21 
December  2006  made  near  the  Athens  DPS-4  digisonde  station  (38.0°  N;  23.5°  E).  In  Figure  1, 
the  radial  projections,  onto  the  Earth’s  surface,  of  the  RO  tangent  points  are  plotted  in  blue;  the 
station  location  is  indicated  by  the  red  star.  For  this  observation,  the  COSMIC  RO 
measurements  were  made  close  to  the  digisonde  location,  with  a  moderately  large  RO  smear 
parameter  of  about  1010  km. 
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Figure  1.  A  COSMIC  satellite  pass  over  the  Athens  digisonde  on  December  21, 2006.  The 
blue  symbols  show  radial  projections  of  the  ray  path  tangent  points.  The  red  star  indicates 
the  station  location. 

Figure  2  presents  the  ionograms  recorded  at  0800  UT  (10  LT)  and  0815  UT.  It  usually  takes 
about  three  minutes  to  complete  an  ionogram  scan  so,  more  precisely,  the  ionograms  presented 
cover  the  periods  0800-0803  UT  and  0815-0818  UT,  respectively.  These  are  very  clean  early 
morning  ionograms  with  unambiguous  F  traces  and  well-defined  F  layer  critical  frequencies, 
foF2.  No  E  layer  trace  is  visible  because  of  poor  signal-to-noise  conditions  in  the  lower 
frequency  range  at  the  Athen’s  site.  These  plots  also  show  the  derived  electron  density  profiles 
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(black  curves),  i.e.,  the  plasma  frequency  as  a  function  of  true  height  calculated  from  the 
ionogram.  In  cases  when  no  E-layer  traces  were  visible,  the  E  region  at  about  100  km  was 
modeled.  The  topside  F  layer  profiles  were  approximated  by  an  a  -Chapman  function  [Reinisch 
and  Huang,  2001]. 


Figure  2.  Athens  ionograms  recorded  during  the  occultation  on  December  21, 2006,  at 
0800  UT  and  0815  UT. 

Figure  3  compares  the  digisonde  profiles  for  Athens  with  the  COSMIC  RO  profiles.  The  two 
ionograms  used  for  this  comparison  bracketed  the  time  of  the  RO  observation,  and  the  measured 
foF2  value  increased  only  slightly  from  5.55  MHz  at  0800  UT  to  5.75  MHz  at  0815  UT.  The 
plots  in  Figure  3  show  the  vertical  profiles  of  electron  density  (top)  and  plasma  frequency 
(bottom). 
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Occultation:  0758-0815  UTC 
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Figure  3.  Comparison  of  the  digisonde  and  COSMIC  RO  profiles  for  Athens.  The  plots 
show  the  electron  density  (top  figure)  and  plasma  frequency  (bottom)  as  a  function  of 
height. 


Overall,  the  RO  profile  is  in  reasonable  agreement  with  the  digisonde  profile,  although  the 
RO  peak  density  is  somewhat  smaller  than  the  measured  digisonde  values.  The  RO  measurement 
of  foF2  was  -0.25  MHz  smaller  than  the  digisonde  measurement  made  at  0800  UT  and 
-0.45  MHz  smaller  than  digisonde  foF2  measurement  at  0815  UT.  The  lower  part  of  the  RO  F2 
profile  is  in  very  good  agreement  with  the  F2  digisonde  profile.  The  digisonde  E  region  profile  is 
modeled  (essentially  using  the  IRI  model),  so  no  comparison  should  be  made.  Directly  above  the 
F2  peak,  the  RO  profile  is  slightly  thinner  than  the  digisonde  profile,  and  the  opposite  is  true 
above  350  km.  It  has  been  shown  that  the  topside  a  -Chapman  profile  is  accurate  only  for  the 
first  -200  km  above  hmF2  [Reinisch  et  al.,  2007;  McNamara  et  al.,  2007],  so  comparisons  above 
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this  height  are  not  meaningful.  For  this  particular  measurement,  the  RO  profile  shows  a 
reasonably  close  agreement  with  the  measured  Ne  profile  over  Athens. 

To  assess  the  RO  performance  for  the  low  latitude  ionosphere,  a  comparison  was  made  with 
the  digisonde  measurements  at  Ascension  Island  (7.9°  S;  345.6°  E).  Figure  4  shows  the  RO  and 
digisonde  profiles  for  December  20,  2006,  at  2330  UT  (2230  LT)  on  the  left  side,  and  the 
measured  nighttime  ionogram  on  the  right.  The  range  and  frequency  spread  in  the  ionogram, 
typical  for  nighttime  at  this  equatorial  anomaly  station,  indicates  small-  and  medium-scale 
irregularities  in  the  F  region.  However,  the  local  foF2  value  was  reasonably  accurately  measured 
as  8.5±0.3  MHz.  The  RO  tangent  point  projection  passed  directly  over  the  station  and  the  smear 
length  was  about  1000  km. 


December  20,  2006 


Occupation:  2319-2340  UTC 

ionogram:  2330 UTC  Frequency  [MHz]  Pi 

Figure  4.  Comparison  of  the  Ascension  Island  profiles  derived  from  digisonde  and 
COSMIC  RO  measurements  (plots  on  the  left)  obtained  on  December  20, 2006.  On  the 
right  is  the  measured  ionogram  demonstrating  spread  F  conditions. 

In  this  case,  the  RO  technique  significantly  underestimates  the  peak  density  by  ~2  MHz  or 
30%  in  terms  of  plasma  frequency,  or  -60%  in  terms  of  Ne.  Clearly,  in  this  situation,  the  RO 
profile  is  not  an  accurate  description  of  the  profile  at  Ascension  Island. 

The  next  example  compares  the  RO  measurements  made  near  another  equatorial  digisonde 
station,  Kwajalein  Island  (9.4°  N;  167.4°  E)  on  December  21, 2006,  at  1220  UT  (2320  LT).  In 
this  case,  the  RO  peak  density  is  very  close  to  the  digisonde  value  (only  0.25  MHz  lower),  but 
the  RO  bottomside  profile  is  completely  wrong.  The  ionogram  shows  a  local  sporadic  E  layer 
with  foEs  =  2.7  MHz,  and  probably  the  presence  of  the  Es  layer  has  affected  the  RO  profile 
inversion. 
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Figure  5.  Comparison  of  the  Kwajalein  Island  profiles  obtained  from  digisonde  and 
COSMIC  RO  measurements  on  December  21,  2006. 


The  last  example  is  for  high  latitude  observations  near  two  Siberian  stations:  Zhigansk 
(66.8°  N;  123.4°  E)  and  Yakutsk  (62.0°  N;  129.6°  E).  Figure  6  shows  the  location  of  the  two 
stations  and  the  tangent  point  footprints  for  an  RO  measurement  on  December  21,  2006.  The 
distance  between  the  digisonde  stations  is  610  km.  The  occultation  of  December  21,  2006,  was 
closer  to  the  Yakutsk  station  (see  Figure  6),  and  the  horizontal  smear  of  this  occultation  was 
1417  km.  Figure  6  illustrates  the  geometry  of  these  measurements,  and  Figure  7  compares  the 
electron  density  profiles  at  the  two  stations  with  the  RO  profile  on  December  21, 2006,  at 
0150  UT  (0950  LT).  The  right  panel  shows  the  ionograms  made  at  the  two  stations  at  this  time. 
Clearly,  the  digisonde  measured  peak  electron  density  is  significantly  lower  at  Zhigansk  than  at 
Yakutsk.  It  is  possible  that  the  Zhigansk  station  was  located  inside  the  mid-latitude  trough  at  that 
time.  The  RO  profile  is  more  similar  to  the  Yakutsk  digisonde  profile  than  the  Zhigansk  profile. 
For  this  particular  measurement,  the  RO  profile  footprint  essentially  extended  longitudinally  (see 
Figure  6),  and  with  the  trough  located  north  of  Yakutsk,  it  probably  did  not  seriously  affect  the 
RO  inversion.  However,  in  situations  with  significant  latitudinal  smear  over  ~1000  km,  the  mid¬ 
latitude  trough  would  likely  complicate  the  RO  inversion.  In  the  present  example,  the  relatively 
thick  E  layer  with  foE  =  2.2  MHz  at  ~150-km  altitude,  as  indicated  by  the  RO  profile,  was  not 
confirmed  by  the  Yakutsk  ionogram. 


7 


December  21  occultation 


GEO  Longitude 


Figure  6.  RO  tangent  points  (blue)  near  Zhigansk  and  Yakutsk  digisonde  stations  on 
December  21,  2006 


December  21,  2006 


lonograms:  0150UTC 


Figure  7.  RO  measurements  on  December  21,  2006,  in  the  vicinity  of  Zhigansk  and 
Yakutsk  digisonde  stations. 
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Finally,  for  all  25  overpass  events  analyzed  (20  and  21  December  2006),  the  foF2  values 
determined  from  the  digisonde  ionograms  and  the  corresponding  RO  measurements  were 
compared. 


Digisonde  foF2  [MHz] 

Figure  8.  Comparison  of  foF2  values  deduced  from  simultaneous  digisonde  and  RO 
measurements. 

Statistically,  the  RO-measured  foF2  values  lie  reasonably  close  to  those  obtained  from  the 
digisonde  observations.  Out  of  25  cases,  there  are  five  foF2  measurements  with  an  error  over 
1  MHz;  all  other  measurements  are  within  this  error  margin.  Larger  errors  appear  to  be 
associated  with  larger  foF2  values. 

2.4  Conclusions  and  Future  Efforts 

The  results  of  these  comparisons  are  interesting,  and  we  now  plan  to  continue  analyzing  data 
for  other  time  periods  and  locations.  For  the  amount  of  data  analyzed,  the  RO  measurements  of 
the  foF2  values  on  average  were  close  to  those  measured  by  collocated  digisondes  with  a  typical 
difference  of  less  than  1  MHz.  There  were,  however,  a  few  cases  in  which  significantly  larger 
differences  were  observed.  These  more  extreme  measurements  appear  to  be  associated  with  the 
equatorial  region.  It  also  should  be  pointed  out  that  for  the  December  21-22,  2006,  campaign,  we 
were  unable  to  make  a  data  comparison  for  the  high-latitude  region  because  the  predicted 
occultations  over  the  digisonde  stations  in  this  region  did  not  produce  ionospheric  profiles.  In  our 
future  work,  we  will  select  several  of  the  most  representative  digisonde  sites  for  which  more 
comprehensive  comparison  will  be  made.  In  the  last  quarter,  we  had  developed  a  software  code 
for  extracting  and  processing  files  in  the  raw  CDF  (common  data  format)  that  are  stored  in  the 
COSMIC  database.  This  will  significantly  speed  up  data  analysis  since  we  will  not  have  to 
manually  select  the  files  from  the  database  relying  on  calculated  satellite  positions  as  in  the  past. 
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With  this  tool,  it  will  be  possible  to  accumulate  enough  data  for  the  comparison  in  equatorial, 
middle,  and  high  latitude  regions.  Our  ultimate  objective  is  to  construct  a  statistical  picture  of  the 
validity  of  the  RO  technique  and  to  establish  its  advantages  and  limitations. 


3.  ARTIST  UNCERTAINTY  SPECIFICATION 

3.1  Introduction 

The  ARTIST  ionogram  autoscaling  software  is  an  intelligent  system  that  UMLCAR 
developed  for  the  extraction  of  ionospheric  specification  data  from  digisonde  ionograms.  There 
is  a  long-standing  need  to  enhance  the  ARTIST  algorithm  in  order  for  the  assimilative 
ionospheric  models  such  as  the  GAIM  [Schunk  et  al.,  2004]  to  know  the  uncertainty  (i.e.,  error 
bars)  associated  with  the  ARTIST-derived  true  height  7V(/i)-profiles.  Previous  efforts  at 
UMLCAR  were  directed  towards  the  development  of  a  sensible  technique  for  automatic 
calculations  of  the  profile  uncertainty  boundaries.  In  the  current  project,  the  developed 
methodology  is  being  implemented  in  the  operating  digisonde  software.  We  also  concentrated  on 
development  of  novel  pre-assimilation  solutions  that  would  help  rule  out  ARTIST  results  that  are 
deemed  to  be  of  low  quality. 

3.2  Approach 

To  determine  the  uncertainty  of  the  N(h)  profiles,  we  constructed  “boundary”  profiles  giving 
the  inner  and  outer  limits  for  each  N(h)  point  calculated  by  the  NHPC  inversion  algorithm 
[Reinisch  and  Huang,  2001].  Figure  9  illustrates  this  approach. 


Figure  9.  Point-by-point  uncertainty  of  the  ARTIST  electron  density  profile  is  obtained  by 
calculating  inner  and  outer  boundaries  enclosing  the  profile.  Both  boundaries  rest  on  five 
anchor  points  whose  uncertainties  are  known  from  statistical  analysis  of  manually 
evaluated  differences  between  automatically  and  manually  scaled  values. 
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The  inner  and  outer  boundary  profiles  are  determined  using  a  set  of  anchor  points  whose 
uncertainties  are  known  from  statistical  analysis  of  prior  data.  Figure  9  shows  the  profile  anchor 
points  hsE  (starting  height  of  E  layer),  foE  (plasma  frequency  at  the  E  layer  peak),  fsFl  (plasma 
frequency  at  the  start  of  the  FI  layer),  foFl  (plasma  frequency  at  the  FI  layer  peak),  and  foF2 
(plasma  frequency  at  the  F2  layer  peak). 

The  anchor  point  uncertainties  vary  with  time,  location,  state  of  the  ionosphere,  and 
operational  state  and  configuration  of  the  sounder.  It  is,  however,  possible  to  capture  the  essence 
of  measurement  uncertainties  through  statistical  analysis.  Therefore,  we  have  been  concentrating 
our  efforts  on  the  statistical  representation  of  ARTIST  errors  by  comparing  ARTIST  results  to 
manually  interpreted  “ground  truth”  data. 

3.3  Statistical  Analysis  of  ARTIST  Performance  Against  Manual  Data 

3.3.1  General  Approach 

We  expected  that,  statistically,  ARTIST  uncertainty  would  be  different  for  different 
digisonde  locations,  geophysical  conditions  (such  as  severity  of  ionospheric  disturbance),  solar 
cycle,  as  well  as  with  time  of  day  and  season.  Thus,  we  envision  a  table  of  typical  uncertainties 
for  each  digisonde  that  ARTIST  software  would  use  to  find  the  appropriate  values,  depending  on 
certain  classification  criteria,  and  then  apply  them  to  describe  the  derived  N(h)-profiles.  The 
lookup  criteria  have  to  be  determined  automatically  by  the  ARTIST  itself.  Obvious  candidates 
that  were  considered  for  classification  of  the  ARTIST  results  into  different  uncertainty  groups 
were  the  ARTIST  version,  digisonde  location  and  model,  hour  of  day,  season,  and  severity  of 
spread  F.  Of  all  these  criteria,  the  spread  F  severity  level  is  determined  using  an  intelligent 
processing  algorithm,  while  other  classification  criteria  such  as  time  of  day  and  season  are  stated 
as  measurement  attributes. 

3.3.2  Spread  F  Detection  Algorithms 

Figure  10  shows  an  ionogram  recorded  at  Jicamarca  on  November  22,  1998  05:00UT  during 
severe  spread  F  conditions  (left);  the  right  panel  of  10  shows  the  “edge  element”  (edgel)  pattern 
calculated  using  the  assumption  that  every  positive  amplitude  gradient  in  the  mixture  of 
overlapping  echoes  is  a  true  leading  edge  of  an  echo. 


11 


JICAMARCA,  JI91J 


1998.11.22(326)05:00:58 


Frequency  [MHz] 


Figure  10.  A  digisonde  ionogram  recorded  during  severe  spread  F  conditions  in  Jicamarca, 
Peru  (left).  If  multiple  echoes  are  resolved  and  presented  as  individual  edgels  (right),  the 
trace  extraction  algorithm  is  overwhelmed  with  numerous  possibilities  of  grouping  edgels 
to  traces. 


In  the  ARTIST-5  software,  the  severity  of  spread  F  is  evaluated  as  part  of  an  ionogram  pre¬ 
analysis  by  simply  counting  the  number  of  edgels  per  frequency  in  the  F  region  (see  Figure  10). 
We  use  an  M-criterion,  where  M  is  a  median  of  {Mi}  and  M\  is  number  of  edgels  found  on 
frequency  i,  M\  >  0.  Table  1  provides  ranges  of  the  M-criterion  used  in  ARTIST-5  to  identify 
four  categories  of  the  spread  F  conditions.  ARTIST  4  and  4.5  do  not  have  the  M-criterion,  but 
they  autoscale  two  standard  spread  F  characteristics,  FF  (frequency  spread)  and  QF  (range 
spread).  Ranges  of  FF  and  QF  for  the  same  four  categories  are  also  given  in  Table  1 . 


Table  1 .  M-Criterion  Values  for  Four  Categories  of  the  Spread  F  Condition 


Quiet 

Moderate 

Severe 

Very  Severe 

ARTIST-5: 

M-criterion 

<4 

Between  4  and 

6 

Between  7  and 

11 

>11 

ARTIST- 

FF  <  0. 1  MHz 

FF  <  0.6  MHz 

FF  <  2.0  MHz 

FF  >  2.0  MHz 

4,4.5,5:  FF  and 
QF  values 

QF  <  5  km 

QF  <  20  km 

QF  <  70  km 

QF  >  70  km 

3.3.3  ARTIST  Confidence  Level  of  Autoscaling  and  Uncertainty 

We  are  still  not  sure  whether  the  ARTIST’S  Confidence  Level  of  Autoscaling  (CLA),  an 
intelligent  merit  check  of  the  ARTIST  processing  results,  can  be  used  to  alter  typical  uncertainty 
bounds  before  they  are  applied  to  the  output  data.  Our  initial  impression  is  that  the  CLA  utility  to 
spot  obvious  autoscaling  blunders  is  well  established,  and  that  this  capability  should  be  used  to 
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rule  out  those  autoscaled  records  completely  instead  of  letting  them  enter  the  assimilation 
process  even  using  a  high  uncertainty.  The  ability  of  CLA  to  reliably  reflect  minor  mistakes  of 
the  autoscaling  analysis  is  yet  to  be  established. 

3.3.4  Calculating  Uncertainty  Bounds 


Two  major  methods  to  deduce  ARTIST  uncertainties  were  used  in  the  course  of  our  study.  A 
quick  estimate  of  the  overall  uncertainty  can  be  obtained  using  a  cumulative  difference  plot,  as 
shown  in  Figure  1 1 . 

Absolute  difference  between  ARTIST  and  manual  foF2 


Grahamstown,  RSA,  2337  ionograms  in  2005-2006 


Figure  11.  Cumulative  difference  function  of  the  absolute  |autoscaled-manual|  differences 
for  a  rapid  evaluation  of  ARTIST  uncertainty.  ARTIST-5  reports  of  foF2  for 
Grahamstown  have  0.15-MHz  uncertainty  at  95%  level. 

The  cumulative  difference  plots  show  the  percent  of  ionograms  with  absolute  errors  less  or 
equal  to  the  abscissa  value;  they  approximate  the  cumulative  probability  distribution  function. 
The  uncertainty  is  determined  at  a  particular  percentile;  the  example  given  in  Figure  1 1  shows 
that  95%  of  all  ionograms  have  fo F2  errors  of  0.1 5  MHz  or  less. 

A  more  informative  uncertainty  estimate  uses  frequency  histograms  that  treat  positive  and 
negative  errors  separately,  as  in  Figure  12. 
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ARTIST  5  Confident  foF2  scaling 


Figure  12.  Uncertainty  bounds  for  foF2  derived  from  the  frequency  histogram  of 
(ARTIST-manual)  differences  for  Boulder  Digisonde  256  ionograms.  For  12,612  ionograms 
scaled  with  confidence  (87%  of  total  14,501  ionograms),  the  uncertainty  bounds  at  5% 
percentile  are  -0.25  and  +0.35  MHz. 

The  uncertainty  bounds  were  derived  from  the  frequency  histograms  of  (ARTIST-manual) 
differences  at  intersections  of  the  histogram  plot  with  the  line  that  divides  all  ditribution  to  95% 
and  5%  of  the  total  ionogram  count  (5%  percentile  uncertainty).  In  the  Figure  12  example,  the 
uncertainty  bounds  are  derived  as  -0.25  and  +0.35  MHz. 

3.3.5  Boulder  Digisonde  Dataset 

We  began  analysis  of  the  A^(A)-profile  anchor  point  uncertainties  with  the  dataset  from  the 
Boulder  digisonde  for  which  we  have  accumulated  almost  two  years  of  manually  interpreted 
ionograms.  Figure  13  presents  manually  scaled  foF2  values  for  the  Boulder  site  during  2004- 
2006.  Anchor  point  uncertainties  were  evaluated  from  this  massive  amount  of  data.  The 
confidence  level  reflects  the  complexity  of  the  ionogram  (ionospheric  tilt,  blanketing  Es, 
absorption,  etc.)  as  well  as  the  ionogram’s  quality. 
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Figure  13.  Availability  of  manually  interpreted  data  for  the  Boulder  Digisonde  station.  The 
plot  shows  foF2  for  every  processed  ionogram  for  years  2004-2006. 

As  a  first  step  in  our  analysis,  the  Boulder  ionogram  measurements  were  divided  into  two 
groups:  one  for  the  winter  period  and  one  for  the  summer  period.  Then,  within  each  period,  the 
uncertainties  for  nighttime  and  daytime  ionograms  were  calculated.  We  then  deduced  the 
ARTIST  uncertainty  using  the  cumulative  difference  characteristic  representing  the  difference 
between  the  ARTIST  and  manually  scaled  foF2  parameter  for  the  same  ionogram.  Figure  14 
shows  the  results  of  the  uncertainties  calculation  for  46,226  Boulder  digisonde  ionograms  scaled 
by  ARTIST  4.5.  Plots  show  foF2  cumulative  differences  calculated  for  nighttime  and  daytime 
ionograms  for  each  period.  The  initial  results  show  that  ARTIST  foF2  uncertainty  (determined 
at  90%  cumulative  difference  level)  is  equal  to  0.31  MHz  during  the  wintertime  measurements 
with  insignificant  difference  between  nighttime  and  daytime  measurements,  and  0.17  MHz  for 
summertime  measurements,  again  with  a  small  difference  between  nighttime  and  daytime 
ionograms. 

3.3.6  Classification  by  Spread  F  Severity  Level  and  CLA 

A  comparison  study  is  underway  to  confirm  the  feasibility  of  detecting  the  level  of  spread  F 
in  ionograms  by  ARTIST-5  in  order  to  apply  different  sets  of  uncertainty  bounds  to  reported 
data.  Figure  1 5  presents  an  example  of  such  a  statistical  study  on  manually  scaled  Gakona,  AK, 
DPS-4  digisonde  data.  It  is  clear  from  the  figure  that,  indeed,  automatically  calculated  level  of 
spread  F  for  each  ionogram  can  be  efficiently  used  as  the  criterion  for  associating  typical 
uncertainty  bounds  on  the  ionogram-derived  data.  In  this  example,  comparing  ARTIST  5.0.2-b7 
scaled  foF2  values  to  manually  interpreted  values,  the  95%  uncertainty  bound  (i.e.,  the  foF2 
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autoscaling  error  value  that  only  5%  of  all  processed  data  exceed)  is  0.33  MHz  for  quiet 
conditions,  0.39  MHz  for  moderate  spread  conditions,  and  0.6  MHz  for  severe  spread  conditions. 


Boulder  Winter  Boulder  Summer 


AfoF2,  MHz 


Figure  14.  ARTIST  4.5  uncertainty  calculated  for  the  winter  and  summer  time  for  Boulder 
digisonde.  Plots  show  foF2  cumulative  differences  calculated  for  nighttime  (black)  and 
daytime  (red)  ionograms  for  each  period. 


Use  of  Confidence  Level  for  foF2  reports 

Gakona  DPS-4,  2002-2006.  Total:  5474,  foF2  detected  in  4569,  confident  data  in  2348  (51%) 
Quiet:  confident  1615  of  2405  (67%),  moderate:  550  of  1513  (36%),  severe:  183  of  651 

(28%) 


Error,  M  Hz 

Figure  15.  Uncertainty  of  foF2  autoscaled  value  for  Gakona,  AK,  digisonde  as  a  function  of 
automatically  detected  level  of  spread  F  conditions.  Comparison  of  ARTIST  5.0.2-b7 
results  for  5474  manually  scaled  ionograms  during  the  2002-2006  period. 
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Thus,  instead  of  applying  universally  0.39-MHz  uncertainty  of  foF2  on  all  ionograms, 
differentiation  by  severity  of  ionospheric  conditions  allows  a  more  correct  report  of  0.33  MHz 
uncertainty  of  foF2  for  ionograms  collected  during  quiet  conditions,  and  0.6  MHz  uncertainty  on 
ionograms  obtained  during  severe  spread  F. 

If  ionograms  with  low  confidence  are  removed  from  the  analysis,  the  uncertainty  bounds 
become  smaller.  While  ~  92%  of  all  ionograms  are  scaled  by  ARTIST-5  confidently  (CLA  > 
50%)  at  Grahamstown,  interpretation  at  the  high  latitude  station  at  Gakona,  Alaska,  is  confident 
only  on  half  of  all  ionograms,  with  only  28%  of  the  ionograms  processed  confidently  during 
severe  spread  F  conditions. 

3.3.7  Uncertainty  Associated  with  Different  Versions  of  ARTIST 

Figure  16  presents  comparison  of  the  ionogram  processing  results  obtained  with  ARTIST4.5 
and  ARTIST5.0  versions  of  the  program. 


AfbF2  =0.27  MHz 
AfoFI  =0.29  MHz 
AfbE  =  0.48  MHz 


AfoF2=  0.30  MHz 
AfoFI  =  0.35  MHz 
AfoE  =  0.55  MHz 


Figure  16.  Cumulative  difference  characteristics  for  Boulder  2004  ionograms.  This  figure 
compares  the  results  obtained  with  ARTIST5.0  left)  and  ARTIST4.5  (right).  Numerical 
values  for  foF2,  foFl,  and  foE  scaling  uncertainties  (calculated  at  95%  level)  are  listed  as 
well. 


The  Boulder  ionograms  for  the  year  2004  were  used  as  input  data  for  this  comparison.  As 
before,  the  accuracy  was  determined  as  95%  uncertainty,  e.g.,  AfoF2  =  0.3  MHz  means  that  for 
95%  of  the  automatically  processed  ionograms,  foF2  values  calculated  by  ARTIST  differ  from 
the  manually  interpreted  ones  by  less  that  0.3  MHz.  In  the  calculation  of  these  characteristics  no 
spread  F  condition  detector  was  applied  since  ARTIST4.5  does  not  have  that  capability.  Clearly, 
the  latest  ARTIST5.0  version  provides  noticeable  improvement  in  the  scaling  accuracy  for  all 
three  ionogram  parameters.  It  is  also  worth  noting  that  for  Boulder  the  accuracy  of  foE  scaling 
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was  significantly  lower  that  that  for  foFl  and  foF2.  This  can  be  attributed  to  a  higher  background 
noise  level  in  the  lower  frequency  range  (<  3  MHz)  characteristic  for  that  station. 

3.3.8  Station-Specific  Uncertainty 

While  we  concluded  that  there  is  little  uncertainty  difference  between  day  and  night  time 
ionograms,  evaluations  for  different  digisonde  locations  show  quite  remarkable  station-specific 
differences  in  the  ARTIST  processing  outcome.  We  performed  statistical  evaluation  of  the 
automatic  scaling  over  a  large  amount  of  data  collected  in  2003-2004  at  the  following  digisonde 
stations:  Boulder,  Ebro,  Athens,  Jicamarca,  and  Grahamstown.  For  all  stations,  except  Boulder, 
the  4.5  version  of  the  ARTIST  software  was  used,  because  it  is  the  version  currently  installed  on 
those  sounders.  Results  of  the  analysis  are  presented  in  Table  2  that  the  Boulder  and  Ebro 
stations  operate  the  older  sounder,  the  DGS256.  These  two  stations  also  offer  significantly  larger 
amounts  of  manually  scaled  data  available  for  the  purpose  of  such  comparison  (see  columns 
titled  “Cases”). 


Table  2.  ARTIST  Uncertainty  for  Various  Digisonde  Locations 


Digisonde 

AfoF2, 

MHz 

Cases 

AfoFl, 

MHz 

Cases 

AfoE, 

MHz 

Cases 

Boulder,  A4.5 

DGS  256 

0.30 

28,609 

0.35 

5,499 

0.55] 

9,999 

Boulder,  A5.0 

DGS  256 

0.27 

10,798 

0.29 

2,029 

0.48 

4,103 

Ebro,  A4.5 

DGS  256 

0.75 

30,746 

0.35 

5,587 

0.35 

11,580 

Athens,  A4.5 

DPS-4 

0.40 

6,451 

0.25 

455 

0.25 

1,124 

Jicamarca,  A4.5 

DPS-4 

0.50 

6,003 

0.50 

948 

0.35 

1,999 

Grahamstown,  A4.5 

DPS-4 

0.30 

2,022 

0.20 

348 

0.20 

544 

The  numbers  in  the  table  suggest  that  for  the  stations  analyzed,  the  typical  accuracy  of  foF2 
and  foFl  scaling  is  of  the  order  of  0.2-0.5  MHz  with  the  exception  of  the  Ebro  station.  This  is  a 
very  impressive  result  demonstrating  the  power  of  the  ARTIST  algorithm.  For  the  Ebro  station, 
for  which  AfoFl  is  equal  to  0.35  MHz  and  AfoF2  is  equal  to  0.75  MHz,  the  problem  is  the 
presence  of  a  known  interferer  in  the  frequency  band  6-7  MHz  which  often  blankets  a  significant 
portion  of  the  ionogram  F2  trace  near  its  cusp  making  automatic  scaling  extremely  difficult. 
Therefore,  the  observed  degraded  ARTIST  performance  was  not  surprising.  Note  that  the 
interferer’ s  presence  does  not  affect  ARTIST  scaling  of  the  FI  trace  which  is  found  at  a  lower 
frequency  range.  The  Jicamarca  DPS  station  presents  another  significant  challenge  for  an 
automatic  processing  because  of  the  irregular  nature  of  the  equatorial  ionosphere.  The  resulting 
uncertainty  of  0.5  MHz  reflects  the  complicated  ionogram  pattern  found  there.  In  the  next  step 
we  will  investigate  the  effect  of  applying  the  spread  F  condition  detector  algorithm  to  the 
Jicamarca  ionogram  processing.  The  accuracy  of  foE  scaling  is  between  0.2  MHz  and  0.5  MHz, 
which  is  a  large  error  given  the  significantly  smaller  absolute  value  of  foE.  We  are  currently 
exploring  several  ways  for  improving  ARTIST’S  performance  in  the  E  region. 
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Our  main  results  of  the  ARTIST  accuracy  analysis  accumulated  thus  far  can  be  summarized 
as  follows: 

•  Statistically,  ARTIST  does  an  excellent  job  (see  Table  2). 

•  There  is  no  pronounced  seasonal  dependence  of  ARTIST  precision. 

•  Station-to-station  variability  of  the  ARTIST  accuracy  is  significant  because  of  the 
differences  in  the  noise  (interference)  environment  at  the  stations  and  levels  of 
complexity  of  the  ionogram  trace  signatures. 

3.4  Reporting  Uncertainties  in  SAO  Format 

A  major  effort  is  underway  to  report  uncertainty  bounds  in  the  output  files  of  the  ARTIST5 
autoscaler.  We  have  been  successfully  using  new  SAOXML  5.0  format  for  the  ionogram-derived 
data  [http://umlcar.uml.edu/SAOXML/]  for  over  two  years.  So  far,  SAOXML  5.0  is  used  to 
report  the  same  amount  of  data  as  in  previous  version  of  SAO  4.3.  The  UMLCAR  software  is 
being  modified  to  read  and  write  uncertainty  information  in  SAOXML  5.0. 

Conceptually,  only  the  ARTIST  software  currently  has  the  capability  to  add  uncertainty 
bounds  to  the  results  of  its  ionogram  analysis.  Modifications  to  SAO  Explorer,  DigisondeLib 
library,  and  the  ionogram  visualization  routines  were  therefore  limited  to  reading/writing, 
internal  organization,  and  proper  display  of  the  uncertainty  data.  The  ARTIST5  software  has 
algorithms  for  calculating  all  necessary  boundaries  from  the  typical  uncertainty  bounds  for  the 
anchor  points  of  the  profile.  As  soon  as  ARTIST5  software  is  modified  to  access  typical 
uncertainty  bounds  individually  for  each  sounder  location  and  level  of  ionospheric  disturbance, 
the  pipeline  of  automated  uncertainty  reports  for  the  digisonde  data  will  be  complete. 

3.4. 1  Uncertainty  Bounds  for  N(h)  Profiles 

Two  styles  of  presenting  uncertainty  bounds  for  the  Ne  electron  density  profile  were 
identified.  First,  each  tabulated  profile  point  can  be  associated  with  the  appropriate  lower  and 
upper  bounds  of  the  reported  electron  density  obtained  from  the  internal  and  external  boundaries 
of  the  profile.  Second,  both  inner  and  outer  boundaries  can  be  stored  analytically  in  terms  of  the 
shifted  Chebyshev  coefficients  representing  the  boundary.  While  we  are  trying  to  establish 
contact  with  the  GAIM  experts  at  Utah  State  University  to  identify  the  preferred  style  of 
reporting  the  boundaries,  both  techniques  are  implemented  in  the  SAOXML  5.0,  at  the  expense 
of  certain  increase  in  the  output  file  size.  The  overhead  of  dual  reporting  of  the  uncertainty 
bounds  for  No  profile  is  considered  at  this  time  insignificant,  considering  the  overall  miniscule 
amount  of  the  ionogram-derived  data  per  ionogram  measurement,  by  contemporary  data  archival 
criteria. 

3.4.2  Software  Modifications 

Modifications  to  the  SAO  Explorer,  DigisondeLib,  and  visualization  libraries  to  support  data 
uncertainties  are  completed  and  tested.  Slight  modifications  to  the  standard  DTD  schema  of 
SAOXML  5.0  were  necessary  to  allow  storage  of  the  uncertainty  bounds  in  the  form  of  profile 
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coefficients.  Figure  17  shows  the  SAO  Explorer  Ionogram  panel  displaying  the  profile 
uncertainty  boundaries  calcualted  for  the  Boulder  ionogram  taken  on  October  6,  2005,  19:45  UT. 

3.5  Future  Efforts 

While  the  SAOXML  5.0  format  has  been  confirmed  and  tested  as  the  vehicle  for  delivering 
the  ARTIST  incertainty  information  to  the  GAIM,  it  remains  to  be  reviewed  and  accepted  by  the 
GAIM  community.  The  review  process  may  result  in  software  modifications. 


Figure  17.  Calculated  Profile  Uncertainty  Boundaries  Shown  in  the  SAO  Explorer. 


The  feasibility  of  using  a  spread  F  severity  detector  to  adaptively  select  uncertainty  bounds 
for  the  reported  ARTIST  characteristics  has  been  proven.  Additionally,  ARTIST5  needs  to  be 
modified  so  as  to  use  an  external  configuration  file  to  obtain  typical  uncertainties  for  a  particular 
location  and  apply  them  to  the  output  data. 

The  task  of  determining  uncertanties  via  statistical  analysis  of  the  autoscaled  versus  manually 
scaled  ionograms  remains  a  labor-intensive  undertaking.  Our  work  on  ARTIST  analysis  will 


20 


continue  into  the  next  year  on  this  project.  We  will  concentrate  on  developing  more  sophisticated 
confidence  level  detectors  that  will  make  it  possible  to  categorize  the  ionograms  based  on  their 
difficulties  of  scaling.  This,  in  turn,  will  allow  deducing  the  associated  measurement  accuracy  for 
the  most  representative  types  of  ionograms. 

4.  DIGISONDE  GROUND-TRUTH  DATA  FOR  THE  NWRA  TOMOGRAPHY 
VALIDATION  PROVJECT 


UMLCAR  continued  to  supply  ground-truth  manually  verified  foF2  and  hmF2  values  to 
support  North  West  Research  Associates,  Inc.  (NWRA)  efforts  to  validate  their  TEC-based 
tomography  calculations  from  GPS  receiver  data.  The  manual  foF2  and  hmF2  values  are  derived 
from  the  digisonde  ionograms  acquired  at  the  Gakona  and  College  sites  during  2001-2006. 

Table  3  describes  the  current  state  of  the  dataset  acquisition  and  processing.  The  total 
number  of  manually  scaled  ionograms  for  this  project  is  15,613. 


Table  3.  Current  State  of  Digisonde  Data  Acquisition  and  Processing  forNWRA  GPS 

Tomography  Project 
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Request  for  digisonde  data  generated  by  NWRA 
+  Dgi sonde  data  acquired  in  □  DBase 
+  Edited  data  reported  to  NWRA 


Total  ionograms  scaled  by  1  September  2007:  1 561 3 


5.  PRECISION  GROUP  HEIGHT  MEASUREMENTS  OF  E  LAYER  VIRTUAL 
HEIGHTS 

5.1  Introduction 

In  the  frame  of  the  present  project,  we  started  work  on  the  development  and  implementation 
of  a  routine  phase  measuring  technique  for  improved  virtual  height  measurements.  Since  the 
development  of  the  first  ionospheric  sounders,  there  have  been  constant  efforts  aimed  at 
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improving  the  accuracy  of  the  measured  characteristics  of  the  ionosphere.  Improved  virtual 
height  accuracy  not  only  improves  the  accuracy  of  the  vertical  electron  density  profiles,  but  is 
also  important  for  monitoring  the  dynamic  processes  in  the  ionospheres:  tides,  gravity  waves, 
sporadic  layers,  neutral  wind  effects,  etc.  It  is  known  that  the  characteristic  vertical  spatial  scales 
of  such  processes  range  from  100  m  to  several  kilometers.  Therefore,  in  order  to  investigate 
these  phenomena,  an  adequate  height  accuracy  is  essential.  We  focused  on  the  virtual  height  of 
the  E  layers  where  the  interaction  between  the  neutral  atmosphere  and  ionosphere  is  the  most 
significant.  Dynamic  coupling  from  the  lower  atmosphere  via  tides  and  waves  plays  an 
important  role  in  the  dynamo  region,  where  energy  is  injected  during  geomagnetic  storms  and 
then  transferred  to  the  higher  heights  in  the  ionosphere. 

For  an  ionosonde,  two  major  parameters  characterize  the  accuracy  of  the  measurements: 
frequency  resolution  and  time  delay  accuracy.  Frequency  resolution  is  simply  determined  by  the 
frequency  step  size  selected  by  the  digital  ionosonde  to  scan  through  the  ionogram  and  in 
principle  can  be  set  to  any  desired  value.  Usually  it  chosen  to  be  20-100  kHz,  which  in  most 
cases  is  sufficient  for  determining  the  features  of  the  ionospheric  electron  density  profile. 
Determining  the  exact  echo  delay  time  is  more  complicated.  As  in  any  radar  system  the 
digisonde  measures  the  time  of  arrival  of  the  echo  pulse  in  relation  to  the  time  of  transmission  of 
the  pulse.  The  pulse  delay  time  r  is  then  converted  to  the  group  or  virtual  height  h'  equal  to 
ct  12  where  c  is  the  free-space  speed  of  light.  The  virtual  height  h'(f)  can  then  be  inverted  to  the 
true  height  using  existing  true  height  inversion  algorithms,  e.g.,  the  Digisonde’ s  NHPC  program 
(Huang  and  Reinisch,  1996).  Obtaining  an  accurate  profile  requires  accurate  h'(f)  values.  The 
usual  ionosonde  technique  of  finding  r  from  the  amplitude  signature  of  the  echo  pulses  has  a 
limited  accuracy  of  several  kilometers,  typically  5-10  km.  The  digisonde’s  phase  technique, 
discussed  below,  which  is  available  in  each  digisonde  but  rarely  applied  because  ARTIST  is  not 
yet  trained  to  make  use  of  the  phase  measurements,  provides  an  unprecedented  accuracy  of  less 
than  1  km  in  measuring  h'(f). 

5.2  The  Virtual  Height  Phase  Technique 

The  technique  for  measuring  the  range  of  the  signal  reflection  using  the  phase  difference 
between  two  separated  frequencies  is  often  referred  to  as  the  d<j>  /  df  method  [Reinisch,  1996; 
Reinisch  et  al.  1997;  Haines,  1994;  Bibl  and  Reinisch,  1978].  Making  use  of  the  phase 
information  carried  by  the  returned  signals  allows  significant  improvement  in  the  accuracy  of 
ionosonde  sounding  in  comparison  to  the  conventional  time  delay  methods  [Davies,  1990].  The 
principles  of  the  d<f>  /  df  method  in  a  free  space  medium  are  straightforward,  but  in  ionospheric 
plasma  the  situation  becomes  more  complicated.  This  section  provides  the  phase  relations  for  the 
signals  propagating  in  an  ionized  plasma  characterized  by  its  refractive  index  n( r,  a>) ,  which  is  a 

function  of  position  r  and  signal  frequency  oo.  For  bottomside  ionospheric  sounding,  the  external 
magnetic  field  B  can  be  considered  constant  with  height.  We  start  with  the  eikonal  solution  for 
plane  waves  in  a  slowly  varying  stratified  medium: 


(1) 
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For  the  case  of  vertical  propagation,  we  can  replace  r  by  h : 

-i(col-[  R  k  dh-tj>^\ 

E(h,t)  =  E0e  [  J*  j=£0e 
where  (f>  =  -[^cot-^*kdh-<f)(^ 

=  [  R  k  dh  +  (f>0  for  t=0, 

J  hfs 


(2) 


where  /io  stands  for  the  height  at  which  the  signal  is  emitted,  e.g.,  ho  =  0  for  groundbased 
sounding.  It  is  important  that  in  transmission,  d<f>0  ldco  =  0 ,  and  digisondes  are  designed  to 
satisfy  this  condition.  Using  the  relation  k  =  0)-n(h,0))/c  (c  is  the  speed  of  light  in  vacuum), 
equation  (2)  is  recast  as: 

rhR  (a))  2  fMa>) 

#*(®)-A=J0  kdh  =  ~l  <o-n(h,co)dh  (3) 

JO  Q  JO 

with  hR  {cd)  denoting  the  height  of  the  signal  reflection.  Differentiating  this  equation  with  respect 
to  the  frequency  we  arrive  at  the  following  relation: 


^-  =  -\ndh=-h\co). 
dco  c '  c 


(4) 


Here,  h'{co)  is  the  virtual  reflection  height  as  a  function  of  sounding  frequency.  That  allows  us  to 
express  the  phase  difference  between  two  frequencies  a)\  and  a>i  as: 


A, J  =  j  =  J  h'(co)dco 

i  dco  J 

<»i  a) 


Using  the  mean  value  theorem,  one  obtains: 

A> -)*'(«>)*»  =  — 

i  c 

with  h'(co^)  <  ht'<  h'(co2) ,  from  where  it  follows  that 


h'= 


4  jdsf 


A  J 


(5) 

(6) 

(7) 


This  expression  makes  it  possible  to  determine  the  reflection  height  h, '  using  the  phase 
difference  A m(f)  obtained  from  the  simultaneous  sounding  at  the  frequencies  co\  and  coi.  There  are 

a  couple  of  complications,  however.  In  actual  ionospheric  sounding  it  is  not  possible  to  directly 
measure  the  total  phase  of  the  signal,  since  the  measured  phase  always  contains  an  unknown 
number  of  multiples  of  2n .  Taking  the  2 n  ambiguity  into  account,  equation  (7)  can  be  written 
as: 
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h'=  —^—kj  =  77 T^Jmeas  +  2^m) 
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or 

c  cA  <b 
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(8) 


(9) 


2A/  4 nAf 

where  A^^  is  the  actual  measured  phase  difference  (which  is  always  in  the  interval  0  <-»  2n), 

and  m  is  an  unknown  integer  multiplier.  This  2n  ambiguity  in  the  phase  measurements  leads  to 
the  ambiguity  in  determining  the  range  of  the  reflection  from  phase  measurements  alone.  It  is 
easy  to  see  from  (9)  that  the  ambiguity  A K  in  the  range  measurement  is 

<l0> 

In  the  digisonde,  the  frequency  separation  A f  is  usually  selected  as  a  multiple  of  5  kHz.  Thus,  if 
A f  =  1- 5 kHz ,  with  /  =  1,2,3...,  then  the  2n  ambiguity  is 


c  3xl08m/j  30  . 

- = - =  —  km 

2Af  21  -5kHz  l 


(11) 


To  resolve  this  ambiguity,  the  group  travel  time  r  determined  from  the  echo  pulse  is  used, 
which  is  always  measured  with  an  accuracy  of  the  order  of  5-10  km.  This  extra  information  is 
used  to  determine  the  first  term  in  (9)  while  the  precise  virtual  height  h'  is  determined  using  the 
sum  of  the  two  terms,  with  the  second  term  providing  the  desired  precision,  which  can  be 
estimated  as 


Sh'  =  8(C^^meas)=^— 
4  nAf  2Af 


)  =  —  £(  A  Jmeas  )  ^ 
2  n  l  2  n 


If  we  assume  that  A  m(j>meas  can  be  measured  with  an  accuracy  of  better  than  2 n  /  32  (which  is 

easily  surpassed  given  for  reasonable  signal-to-noise  ratios),  then  the  h'  measurement  precision 
will  be  at  least 


5 ht' «1/ /  km  (13) 

i.e.,  for  5-kHz  frequency  separation  it  is  better  than  1  km. 


Another  complication  involves  the  difference  in  the  reflection  heights  h'(a>x)  and  h'(a>2)  of 
the  signals  at  the  two  frequencies.  From  equation  (6)  it  is  evident  that  h]  is  bounded  by  h'(co{  ) 
and  h\co2) ,  and  therefore  the  phase  technique  has  a  limited  IT  resolution  resulting  in  an 

uncertainty  u  of  u  <  ^\h' (a>2)  - h' (a>x)\ .  Therefore,  the  frequency  difference  <x>2  -  o\  should  not  be 

too  large.  We  estimated  the  height  differences  using  routine  ionogram  recordings.  For  the  mid- 
latitudinal  ionosphere  (Millstone  Hill)  during  daytime,  the  figures  given  in  Table  4  characterize 
the  uncertainty  of  the  precision  group  height  measurements  for  signals  separated  by  5  kHz. 
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Table  4.  Typical  Precision  Group  Height  Uncertainties 


Region 

Uncertainty  (5  kHz) 

E 

25  m 

E  cusp 

125  m 

F  bottom 

75  m 

F  cusp 

750  m 

The  numbers  in  Table  4  show  that  the  uncertainty  is  large  only  near  the  F  layer  cusp.  Even  at  the 
F  cusp,  however,  the  uncertainty  is  smaller  than  1  km 


5.3  Measuring  E  :Layer  Heights  at  Millstone  Hill 


A  multi-frequency  precision  group  height  (MPGH)  technique  has  been  implemented  into  the 
DPS-4  sounder  operating  at  the  Millstone  Hill  Observatory.  In  September  and  November  2006, 
the  digisonde  was  programmed  to  carry  out  routine  MPGH  measurements  of  E  layer  heights, 
operating  in  the  so-called  “drift  mode”.  In  total,  there  were  40  days  of  observations  accumulated 
with  a  measurement  cadence  of  15  minutes  made  in  the  intervals  of  1  -  21  September  and  3 1 
October  -  19  November  2006.  The  purpose  of  this  campaign  was  to  monitor  the  dynamics  in  the 
E-region  and  to  evaluate  the  performance  and  robustness  of  the  MPGH  technique.  The  following 
four  frequencies  were  used  for  sounding  in  the  frequency-multiplexing  mode:  2330  kHz, 

2335  kHz,  2340  kHz,  and  2345  kHz.  Making  use  of  the  well-known  expression  relating  the 


plasma  frequency  fp  and  the  electron  density  Ne  (  Ne  = 


1 

80.6 


fp  ,  where  Ne  is  in  m‘3  and  /  is 


in  Hz)  these  frequencies  provided  measurements  at  the  heights  corresponding  to  electron 
densities  around  6.74xl010m'3.  Only  echoes  from  E  region  heights  between  90  and  135  km  were 
considered.  New  software  codes  for  automatic  processing  of  these  new  measurements  has  been 
developed. 


The  group  of  four  sounding  frequencies  produced  three  different  frequency  pairs  with  5  kHz 
separation.  For  each  sounding  frequency  we  collected  echoes  from  1 8  ranges  spaced  by  5  km  at 
the  four  antennas.  After  Fourier  transforming  the  received  signals  digisonde  software  selects  the 
strongest  spectral  components  for  processing.  Only  signals  at  least  12  dB  above  the  noise  level  of 
a  spectrogram  were  included  in  the  analysis.  The  complex  Fourier  transform  also  gives  the  signal 

phases  from  which  the  phase  differences  A„,^ _ used  in  (9)  were  calculated.  A  number  of  tests 

were  used  to  improve  the  reliability  of  the  results:  only  signals  that  have  similar  phases  on  all 
four  antennas  were  considered  in  order  to  discard  off-vertical  echoes;  also,  only  those  calculated 
precision  group  heights  were  included  that  were  within  7  km  of  the  h'  value  determined  by  the 
coarse  time  delay  measurement  using  only  echo  amplitudes.  Thus,  for  each  digisonde  sounding, 
it  was  possible  to  calculate  a  large  number  of  independent  precision  group  heights:  one  for  each 
antenna,  each  frequency  pair,  each  coarse  range  and  for  several  of  the  strongest  spectral  lines. 
This  made  it  possible  to  obtain  a  statistical  distribution  of  measured  group  heights  for  a  single 
observation.  Analysis  has  shown  that,  in  most  cases,  a  significant  majority  of  the  calculated 
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group  heights  (>65%)  were  concentrated  within  a  single  statistical  height  bin  1  km  wide  (Figure 
18). 


75n 


Millstone  Hill 
2006.09.14 
14:22  UT 


3  Measurements 
Into  Statistics 


95.5  97.5  99.5  101.5  103.5  105.5  107.5  109.5  111.5  113.5  115.5 

Virtual  height  (km) 

Figure  18.  Distribution  of  virtual  heights  during  one  sounding  made  on  14  September  2006 
at  1422  UT.  A  majority  of  the  calculated  heights  are  concentrated  within  a  single  statistical 
range  bin  at  107.5±0.5  km. 

Given  this  distribution,  it  seems  justified  to  use  the  median  as  the  most  representative  value 
for  the  precision  group  height.  This  processing  algorithm  made  it  possible  to  reduce  a  large 
amount  of  the  raw  sounding  data  to  a  single  precision  group  height  for  each  of  the  observation 
made  once  every  15  minutes. 

Results  of  the  MPGH  measurements  carried  out  in  September  2006  are  presented  in  Figure 
19  showing  the  time  evolution  of  E  layer  virtual  heights.  Naturally,  the  radiowave  reflections 
from  the  regular  E  layer  were  observed  only  during  daytime.  The  presented  data  clearly 
demonstrates  an  expected  daily  pattern  of  the  E  region  virtual  height  variation  with  a  minimum 
height  near  local  noon  (1630  UT  approximately)  and  apparent  maximum  heights  at  sunrise  and 
sunset  (approximately  12  UT  and  21  UT  respectively).  For  these  20  days  of  observations,  the 
average  daily  variation  is  shown  in  Figure  20. 
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Figure  19.  Daily  E  region  precision  group  heights  at  Millstone  Hill  as  function  of  time.  The 
horizontal  axes  represent  time  (Universal  Time),  and  the  vertical  axes  represent  virtual 
height  (km).  Dashed  circles  highlight  significant  decreases  in  E  layer  virtual  heights  and  the 
black  arrows  point  to  hook-shaped  disturbances. 

From  the  monthly  mean  plot,  one  can  see  that  after  sunrise  the  recorded  E  layer  virtual  height 
was  117  km,  which  then  decreases  to  105  km  at  local  noon,  and  subsequently  rises  back  up  to 
122  km  right  before  sunset.  The  calculated  standard  deviations  are  about  4-5  km  near  sunrise 
and  sunset  times,  and  1 .5-2  km  near  noontime.  In  addition  to  the  monthly  average  pattern,  it  is 
possible  to  observe  a  day-to-day  variability  (Figure  20).  Comparing  the  monthly  average  pattern 
in  Figure  20  with  the  daily  variations  recorded  in  September  2006  (Figure  19),  one  notices 
significantly  lower  heights  on  September  2,  5,  8,  and  1 1  (between  14-20  UT),  with  the  measured 
virtual  height  values  as  low  as  98  km.  These  features  were  indicated  in  Figure  19  with  the  dashed 
circles.  The  observed  differences  of  5-7  km  between  the  heights  measured  on  those  days  and 
average  values  are  significantly  larger  than  the  standard  deviations  for  the  corresponding  time 
interval. 
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Time  (UT) 


Figure  20.  Average  daily  variations  of  the  h’E  obtained  from  the  20  days  of  September. 
Vertical  error  bars  indicate  the  standard  deviation  obtained  at  a  given  time. 

Interestingly,  these  ‘low  daily  values’  occurred  with  a  periodicity  of  three  days,  which  could 
be  indicative  of  some  oscillation  activity  at  E-region  altitudes.  Similar  oscillations  have  been 
found  in  the  electron  density  variations  in  the  E  layer  by  other  authors.  For  example,  Zhou 
[  1 998]  reported  a  clear  two-day  modulation  of  the  electron  density  at  the  altitude  range  of  84-94 
km,  which  he  attributed  to  quasi  two-day  planetary  waves.  The  h’E  variations  we  observe  in  this 
work  closely  resembles  height  variations  at  fixed  electron  density  presented  in  Zhou. 

Results  of  the  data  analysis  demonstrate  the  great  potential  offered  by  ionospheric  sounding 
with  enhanced  height  resolution.  It  can  be  concluded  that  in  practice  the  height  precision  offered 
by  the  new  method  is  better  than  1  km.  The  feasibility  and  stability  of  these  measurements  are 
confirmed  by  the  very  narrow  statistical  distribution  of  the  independently  determined  group 
heights  for  each  sounding.  The  accumulated  data  made  it  possible  to  observe  very  interesting 
effects  in  the  daily  variations  of  the  E-region  heights.  We  were  able  to  observe  day-to-day 
variability,  which  in  the  case  of  the  September  measurements  demonstrated  an  interesting  three- 
day  periodicity.  This  dataset  also  shows  repetitive  “hooks”  in  the  h'  records.  It  is  important  to 
point  out  that  such  results  could  not  have  been  obtained  using  the  conventional  time  delay 
measurements  because  the  characteristic  amplitude  of  these  unusual  variations  is  about  one 
kilometer.  Our  preliminary  study  of  the  observed  effects  suggest  that  day-to-day  variability  of 
the  E  layer  height  variations  could  be  related  to  planetary  wave  modulation  of  metallic  ion 
transport  while  the  “hook”  variations  (see  arrows  in  Figure  19)  may  be  caused  by  tidal/gravity 
wave  activity.  The  existing  worldwide  digisonde  network  offers  a  great  advantage  of  organizing 
unprecedented  simultaneous  observations  on  a  global  scale. 
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6.  PROCESSING  OF  DIGISONDE  DRIFT  MEASUREMENTS 


6.1  Introduction 

We  have  been  working  on  the  analysis  of  ionospheric  plasma  drift  velocity  data  accumulated 
in  the  DriftBase  database  over  the  last  few  years.  A  number  of  digisondes  routinely  measure 
ionospheric  “HF  drift”  after  each  ionogram  observation.  It  is  well  known  [Bullett,  1994; 

Reinisch  et  al.,  1998]  that  these  high  frequency  (HF)  drift  measurements  do  not  directly  measure 
the  ion  drift,  but  rather  give  the  changes  in  the  location  of  ionospheric  structures  with  time,  i.e., 
an  apparent  plasma  drift  velocity.  If  these  ionospheric  motions  are  controlled  primarily  by 
plasma  transport  processes  (which  is  usually  the  case  for  middle-  and  high-  latitudes),  then  the 
two  velocities  are  equivalent.  By  comparing  the  digisonde  drift  velocities  to  the  ISR  measured 
ion  drift  at  Millstone  Hill  (a  mid  latitude  station),  Bullett  [1994]  experimentally  confirmed  the 
validity  of  the  digisonde  drift  method.  The  observed  agreement  between  the  two  techniques  was 
best  at  night  and  when  no  traveling  ionospheric  disturbances  (TIDs)  were  present.  Certain 
differences  in  the  two  velocities,  however,  were  observed  during  the  daytime  and  under  non¬ 
equilibrium  conditions  [Reinisch  et  al.,  1998].  In  the  equatorial  region,  this  technique  should 
provide  meaningful  measurements  for  the  horizontal  plasma  velocity  component,  which  is  not 
significantly  affected  by  chemical  processes. 

Currently,  at  most  of  the  digisonde  stations  running  the  drift  mode,  the  plasma  drift  velocities 
are  calculated  automatically  with  the  Digisonde  Drift  Analysis  (DDA)  software  package  and  then 
archived  in  the  UMLCAR  drift  velocity  database  “DriftBase”  (http://ulcar.uml.edu/Drift- 
X.html).  Using  statistical  analysis  of  the  available  drift  data  we  wanted  to  establish  and 
understand  the  typical  daily  “plasma  drift”  variations  as  measured  by  the  digisondes.  We  are 
interested  in  finding  correlations  between  the  behavior  of  the  plasma  velocity  and  other 
geophysical  processes  controlling  the  ionospheric  dynamics. 
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6.2  Data  Analysis 

We  started  a  statistical  analysis  of  the  routine  F  region  digisonde  drift  measurements  using 
data  from  the  Ebro  station  in  Spain  (40.8°N,  0.5°E).  In  total,  we  have  analyzed  18  months  of  the 
vertical  (Vz)  velocity  component  data  collected  during  years  2004  and  2005.  First,  we 
determined  the  typical  daily  pattern  in  terms  of  the  main  harmonics  (diurnal,  semidiurnal,  and 
terdiumal).  The  sum  of  the  first  three  daily  harmonics  reproduces  the  observed  monthly  median 
quite  well  (Figure  21),  except  for  the  deep  negative  excursion  near  sunrise.  The  spectral 
characteristics  of  the  harmonic  components  have  been  computed  from  the  monthly  median 
values  by  a  general  least-square  method. 


Figure  21.  Daily  pattern  of  the  vertical  velocity  data  from  the  Ebro  digisonde  for  December 
2004.  The  15-minute  median  values  of  Vz  are  shown  as  gray  dots.  The  vertical  error  bars 
indicate  the  range  of  80%  of  Vz  values  for  each  given  time.  The  solid  line  represents  the 
daily  pattern  calculated  as  a  sum  of  the  primary  diurnal  harmonics.  The  vertical  dashed 
lines  indicate  the  sunset  (SS)  and  sunrise  (SR)  above  the  station  at  an  altitude  of  300  km. 

We  have  analyzed  the  daily  patterns  of  Vz  and  its  variability  as  a  function  of  local  time  and 
season.  Figure  22  presents  amplitudes  and  phases  of  the  three  main  daily  harmonics  as  a  function 
of  season.  Figure  23  shows  the  seasonal  variation  of  the  phase  of  the  semidiurnal  harmonic  (left 
panel)  and  semidiurnal  harmonic  phase  as  a  function  of  sunrise  time  (right  panel).  The  main 
results  are  summarized  as  follows:  (1)  Vz  displays  larger  variability  during  nighttime  and  during 
winter;  (2)  the  diurnal  harmonic  is  dominant  with  the  larger  amplitudes  observed  at  equinox 
compared  to  that  at  solstice  (Figure  2,  left  panel);  (3)  the  amplitude  of  the  semidiurnal  variation 
is  enhanced  near  equinox,  while  terdiumal  variation  is  strongest  during  the  summer  half-year 
(from  March  to  September)  when  it  has  a  larger  amplitude  than  the  semidiurnal  (Figure  22,  left 
panel);  (4)  the  phases  of  the  diurnal  and  terdiumal  harmonics  are  practically  constant  for  the 
whole  year  (Figure  22,  right  panel);  (5)  the  phase  of  the  semidiurnal  harmonic  is  the  key 
parameter  for  matching  the  minimum  values  of  the  daily  pattern  of  Vz  (Figure  23,  left  panel) 
because  it  is  sun-synchronous,  i.e.,  its  variation  follows  the  change  in  the  sunrise  time  (Figure  23, 
right  panel). 
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Figure  22.  Seasonal  variations  of  the  amplitudes  (left)  and  phases  (right)  of  the  daily 
pattern  of  the  Vz  in  the  F  region.  The  top  plots  show  the  diurnal  component,  the  middle 
and  bottom  ones  correspond  to  the  semidiurnal  and  terdiurnal  component,  respectively. 


Figure  23.  The  left  panel  shows  the  seasonal  variation  of  the  phase  of  the  semidiurnal 
harmonic  (full  line  with  full  dots),  the  time  of  the  minimum  values  of  the  daily  pattern  of 
Vz  (dashed  line  with  open  squares)  and  the  sunrise  time  at  300-km  height  (full  line  with 
open  circles).  The  right  panel  shows  the  scatter  plot  of  the  semidiurnal  phase  against  the 
sunrise  time  for  a  given  month.  Solid  line  is  the  best  linear  fit  with  the  correlation 
coefficient  indicated. 
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The  large  variability  of  Vz  during  nighttime  and  its  winter  enhancement  can  be  in  part 
explained  by  the  neutral  wind  dynamics  and  auroral  energy  inputs.  Other  physical  mechanisms 
may  also  account  for  some  of  the  the  variability.  According  to  the  observed  results,  we  believe 
that  it  is  feasible  to  model  the  daily  pattern  of  the  Vz  component  in  the  F  region  using  digisonde 
drift  measurements.  The  results  presented  here  are  representative  only  for  the  Ebro  station,  and 
they  are  likely  to  vary  with  latitude  and  longitude.  The  formulation  in  terms  of  Fourier  harmonics 
is  quite  simple  (requires  one  parameter  for  each  month)  and  is  easily  updateable.  In  addition,  it  is 
straightforward  to  extend  this  formulation  to  a  larger  geographical  region  by  analyzing  digisonde 
drift  measurements  at  other  latitudes  and  longitudes. 

Larger  data  sets  are  required  to  test/verify  a  possible  solar  activity  dependence,  and  to  assess 
whether  the  proposed  formulation  is  useful  for  global  and  long-term  applications.  Note  that  the 
formulation  in  terms  of  Fourier  harmonics  does  not  precisely  capture  the  Vz  behavior  around 
sunrise.  The  typical  daily  pattern  of  Vz  and  its  seasonal  trend  show  larger  amplitudes  during 
equinox,  agreeing  with  the  expected  seasonal  pattern  as  predicted  by  the  Richmond  et  al.  [1980] 
using  his  empirical  drift  model  (derived  from  ISR  measurements)  This  also  applies  to  the 
variation  in  the  time  of  the  maximum  Vz. 

The  diurnal  component  of  the  Vz  drift  pattern  is  dominant,  calculated  using  the  Ebro 
measurements,  on  the  other  hand,  the  Richmond  model  predicts  a  dominant  semidiurnal 
harmonic.  This,  in  principle,  may  be  a  manifestation  of  the  fact  that  digisonde  measures  an 
“apparent”  plasma  velocity,  while  the  Richmond’s  model  was  developed  for  the  ExB  drift. 
However,  the  daily  Vz  pattern  that  we  have  established  in  this  work  is  an  average  characteristic, 
and  therefore  the  effect  of  irregular  phenomena  like  TIDs  should  not  be  significant.  In  the 
future,  we  will  investigate  this  observed  difference  between  the  model  and  the  measurements  in 
greater  detail.  We  will  concentrate  on  the  role  of  the  neutral  wind  which  has  a  systematic 
behavior  and  may  effect  the  average  vertical  motion  of  the  ionosphere  at  mid-latitudes.  Some 
preliminary  results  of  this  analysis  were  reported  in  an  earlier  report. 

Another  station  for  which  this  statistical  data  analysis  has  been  run  is  the  Jicamarca 
equatorial  station  (12.0°S,  283.2°E).  Figure  24  presents  the  averaged  monthly  patterns  for  the 
digisonde  east-west  velocity  component  at  Jicamarca  in  2006.  The  months  are  indicated  on  each 
panel.  These  data  were  processed  using  the  Drift-Explorer  software  package.  The  monthly  east- 
west  velocity  pattern  is  remarkably  well  defined  and  the  month-to-month  variability  is  also 
apparent.  Error  bars  indicate  standard  deviations  which  are  reasonably  small  for  this  kind  of 
measurements.  Positive  velocity  in  these  plots  corresponds  to  eastward  plasma  flow.  The  change 
from  eastward  to  westward  direction  of  the  plasma  motion,  easily  noticeable  on  the  plots, 
approximately  corresponds  to  the  transition  from  nighttime  to  daytime. 
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Figure  24.  Monthly  averaged  east-west  drift  velocity  data  measured  at  Jicamarca  station 
in  2006. 


Figure  25  presents  the  results  of  spectral  analysis  of  the  east-west  velocity,  namely  the 
monthly  amplitude  spectra.  Not  surprisingly,  for  all  the  months  shown,  the  diurnal  line  is  the 
strongest,  however,  occasionally  the  semidiurnal  and  terdiumal  lines  are  also  appear.  There  also 
appears  a  large  variability  to  the  intensity  of  these  harmonics. 
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Figure  25.  Amplitude  spectra  of  the  east-west  plasma  velocity  measured  at  Jicamarca  in 
2006.  Bottom  axis  shows  the  frequency  in  pHz  while  the  top  axis  shows  the  corresponding 
time  periods  in  hours. 

In  principle,  we  possess  enough  data  for  establishing  a  statistical  model  for  the  east-west 
plasma  drift  velocity  at  Jicamarca  similar  to  the  Richmond  et  al.  [1980]  model,  however,  our 
interest  is  in  analysis  of  the  physical  significance  of  the  observed  variations  and  search  for  a 
possible  correlation  with  other  processes  in  the  ionosphere. 

6.3  Future  Efforts  on  Drift  Data  Analysis 

We  believe  that  these  results  give  evidence  that  the  digisonde  drift  measurements  can  be  used 
for  monitoring  the  F  region  ion  drift.  The  systematic  analysis  of  such  observations  can  contribute 
to  a  better  understanding  of  the  ionospheric  dynamics.  Our  future  efforts  will  extend  this  work  of 
digisonde  F  region  drift  analysis  to  other  digisonde  locations.  We  plan  to  process  more  data  from 
Jicamarca  as  well  as  from  other  stations  using  spectral  analysis  (as  it  was  done  for  Ebro  station 
data)  in  order  to  be  able  to  analyze  seasonal  variations  of  the  plasma  drift  velocity  pattern  and  to 
investigate  its  relationship  with  other  geophysical  processes. 
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6.4  DDA  Software  Upgrade 

To  make  the  drift  data  processing  mor  efficient  we  worked  on  upgrading  the  DDA  software. 
A  translation  of  the  first  half  of  the  existing  Fortran-based  code  to  the  Java  platform  is 
completed.  A  number  of  optimizations  were  made  and  several  programming  errors  were 
corrected  in  the  process.  This  is  a  significant  step  forward  enabling  the  unification  of  the  drift 
data  processing  algorithms  with  the  rest  of  the  digisonde  programs  and  making  use  of  the 
extensive  library  of  the  Java  classes  developed  earlier.  It  is  now  possible  to  easily  maintain 
and/or  update  these  programs,  a  process  that  previously  required  significant  time  consuming 
efforts  because  of  the  difficulty  of  deciphering  Fortran  code.  After  this  translation,  the  Drift- 
Explorer  package  has  an  embedded  capability  of  calculating  the  skymaps  from  the  raw  drift  data, 
which  is  useful  for  both  digisonde  system  testing  and  data  analysis.  The  second  half  of  the  DDA 
software  involving  the  velocity  calculation  routines  will  be  converted  into  Java  in  the  upcoming 
year. 

7.  STUDY  OF  THE  IONOSPHERIC  RESPONSE  TO  STRONG  GEOMAGNETIC 
STORMS 

7.1  Introduction 

The  response  of  the  Earth’s  ionosphere  to  strong  geomagnetic  disturbances  is  still  one  of  the 
primary  interests  in  ionospheric  science  [Rishbeth  et  al.,  1996].  Geomagnetic  disturbances  are 
generated  by  changes  in  the  solar  wind  and  interplanetary  magnetic  field  (IMF)  interacting  with 
the  earth’s  magnetosphere.  The  earth’s  reaction  to  such  perturbations  is  commonly  referred  to  as 
an  ionospheric  storm.  Typically,  these  storms  are  categorized  into  two  characteristic  types: 
positive  storms  (with  the  ionospheric  electron  density  increase)  and  negative  storms  (with  the 
electron  density  decrease).  There  is  fair  agreement  concerning  the  mechanisms  of  the  negative 
storms,  which  are  attributed  to  the  composition  changes  of  the  neutral  atmosphere  resulting  into 
a  faster  recombination  rate,  thus  depleting  the  plasma  [e.g.,  Buonsanto  1999;  Mendillo,  2006].  A 
more  complicated  story  involves  the  positive  storms.  A  number  of  different  mechanisms  have 
been  suggested  to  explain  storm  morphologies  but  a  complete  picture  has  not  yet  emerged. 

Prolss  [1993]  has  offered  the  most  comprehensive  description  of  the  storm  development,  placing 
the  crucial  role  on  the  neutral  wind  disturbances  generated  by  heating  of  atmospheric  gas  at 
auroral  latitudes.  In  his  model  such  neutral  wind  perturbations  are  often  carried  to  lower  latitudes 
by  traveling  atmospheric  disturbances  (TAD)  which  are  a  superposition  of  atmospheric  gravity 
waves.  TADs,  propagating  from  higher  latitudes,  carry  along  meridional  winds  that  push  the 
ionospheric  plasma  upward  along  the  magnetic  field  line  that  in  turn  leads  to  the  increase  in  the 
plasma  density  because  of  the  slower  recombination  rate  at  higher  altitudes  [e.g.,  Rishbeth, 

1998].  This  seems  to  be  a  promising  mechanism  and  its  validity  has  been  confirmed  by  global 
scale  modeling  [Fuller-Rowell  et  al.,  1994].  However,  experimental  verifications  are  still  scarce. 
One  of  the  difficulties  is  that  the  measurements  needed  have  to  be  made  on  a  large  spatial  scale 
and,  in  addition,  measurement  of  the  height  changes  are  also  required.  The  problem  is  that  most 
of  the  ionosphere  measuring  instrumentation  either  does  not  measure  the  height  distribution  of 
electron  density  (e.g.,  GPS  TEC  measurements),  or  cannot  provide  significant  spatial  extent  and 
routine  observations  (e.g.,  incoherent  scatter  radars).  The  large  number  of  existing  ionosondes 
make  it  possible  to  carry  out  this  analysis,  but  many  sounders  are  not  capable  of  automatic 
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calculation  of  the  true  height  distribution  of  electron  density  (i.e.,  inverting  the  measured  virtual 
heights);  therefore,  such  a  task  requires  an  unreasonable  amount  of  time  and  considerable 
human-based  processing.  In  this  situation,  an  existing  worldwide  digisonde  [Reinisch,  1996] 
network  can  provide  the  dataset  suitable  for  such  an  investigation.  It  offers  routine  automatic 
round-the-clock  ionospheric  height  profile  measurements  over  a  wide  range  of  latitudes  and 
longitudes.  Even  though  manual  verification  of  the  ionogram  processing  is  often  recommended, 
the  whole  process  is  greatly  simplified  by  using  the  available  digisonde  software  packages. 

We  took  advantage  of  the  large  amount  of  data  collected  over  a  significant  period  of  time  and 
ran  an  analysis  of  the  mid-latitude  ionospheric  response  to  geomagnetic  storms  using  the 
digisonde  observations  at  different  longitudes  and  latitudes.  For  the  major  storms  of  2001-2005 
we  analyzed  the  evolution  of  the  F  layer  making  use  of  the  electron  density  profiles  and  main 
ionospheric  characteristics  such  as  foF2  and  hmF2.  In  this  reported  research  we  focused  on  the 
very  strong  ionospheric  storm  of  August  24,  2004. 


Table  5.  Digisonde  Station  Coordinates 


Station 

Geo  Coordinates 
(Lat;  Long) 

GM  Coordinates 
(Lat;  Long) 

Europe 

Chilton 

51.5°  N;  359.4°  E 

53.7°  N;  84.4°  E 

Ebro  (Pruhonice) 

50.0°  N;  14.6°  E 

49.6°  N;  98.6°  E 

Roquetes 

40.8°  N;  0.5°  E 

43.2°  N;  81.3°  E 

San  Vito 

40.6°  N;  17.8°  E 

39.9°  N;  98.3°  E 

Athens 

38.0°  N;  23.5°  E 

36.4°  N;  103.0°  E 

America 

Millstone  Hill 

42.6°  N;  288.5°  E 

52.9°  N;  0.3°  E 

Boulder 

40.0°  N;  254.7°  E 

48.3°  N;  320.5°  E 

Wallops  I 

37.9°  N;  284.5°  E 

48.1°  N;  355.6°  E 

Eglin 

30.4°  N;  273.2°  E 

40.3°  N;  342.9°  E 

Analysis  was  limited  to  the  American  and  European  sectors  where  there  exists  a  dense 
network  of  sounders.  Most  of  the  time  the  digisonde  stations  operated  at  a  15-min  cadence. 
Coordinates  of  the  mid-latitude  stations  are  given  in  Table  5.  Unfortunately,  during  strong 
geomagnetic  disturbances  high  latitude  and  polar  stations  often  do  not  produce  useful 
measurements  because  of  strong  radiowave  absorption  in  the  lower  ionosphere. 

7.2  Data  Analysis 

In  Figure  26,  we  see  the  variation  of  the  geomagnetic  parameters  for  the  strong  geomagnetic 
storm  that  occurred  on  August  24,  2005  with  storm  commencement  (SC)  at  0613  UT.  In  this 
figure  the  solar  wind  plasma  density,  plasma  speed,  interplanetary  magnetic  field  (IMF) 
components,  and  SYM-H  index  (DST  index  analog)  are  shown.  The  storm  commencement  at 
0613  UT  coincides  with  a  rapid  increase  in  SYM-H  index  and  similar  abrupt  changes  in  IMF  and 
solar  wind  parameter  values. 
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Figure  26.  Geomagnetic  storm  of  24  August  2005.  The  top  panel  shows  SYM-H  index 
variations,  and  the  next  panels  show  the  solar  wind  plasma  density,  plasma  speed,  and  the 
last  two  panels  show  the  interplanetary  magnetic  field  components.  The  solar  wind  and 
IMF  data  were  taken  using  the  ACE  satellite  and  were  shifted  by  40  minutes  to  account  for 
the  position  of  the  satellite. 

Note  that  at  about  0900  UT  there  is  another  shock  event  accompanied  by  a  southward  turning 
of  IMF  Bz  component.  At  about  the  same  time  the  main  phase  of  the  storm  begins  with  the 
SYM-H  index  reaching  its  minimum  at  1 1  UT.  The  value  of  the  minimum  of  -200  nT  indicates  a 
moderately  strong  geomagnetic  storm 

We  first  compared  the  ionosphere  behavior  on  the  day  of  the  storm  event  to  the  average  quiet 
day  pattern.  Figure  27  introduces  variations  of  the  ionospheric  electron  density  height  profile  as  a 
function  of  time.  In  this  figure  the  left  panels  shows  the  average  quiet  day  data  for  the  European 
stations  Chilton,  Pruhonice,  Ebro,  San  Vito,  and  Athens  with  the  station  locations  given  in  Table 
5.  The  profiles  shown  are  plotted  in  terms  of  the  true  height,  which  is  derived  from  the  original 
measurements  with  the  standard  digisonde  inversion  routine  NHPC.  All  of  the  ionogram  data 
have  been  reviewed  before  being  used  in  the  profile  calculation.  The  topside  parts  of  the  profiles 
are  approximated  by  an  a  -Chapman  function  [Reinisch  and  Huang,  2001].  The  plots  are  in 
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universal  time,  but  ionospheric  sunrise  and  sunset  times  are  indicated  for  one  of  the  stations  with 
red  and  black  semi-circles  respectively,  helping  to  distinguish  daytime  and  nighttime  periods. 


I 

i 

i 

i 


Aug  2005  (Quiet  day  average) 


00  02  04  06  08  10  12  14  16  18  ^  20  22  24 


Universal  Time 


Ai^g  24,  2005  (SSC  0613  UT) 


00  02  04  oq  08  10  12  14  16  18  20  22  24 

Universal  Time 


Figure  27.  Electron  density  profile  variations  as  a  function  of  time  for  the  listed  European 
stations  (see  Table  5).  The  right  panel  shows  storm  day  (August  24,  2005)  data,  while  the 
left  panel  shows  the  quiet  day  average  data.  Electron  density  is  presented  in  terms  of 
equivalent  plasma  frequency  (color  bar).  Sunrise  and  sunset  times  are  indicated  for  the 
Ebro  station  with  red  and  black  circles  correspondingly. 


Data  from  the  different  stations  are  arranged  according  to  the  station’s  geomagnetic  latitude 
as  indicated  on  the  plots.  The  quiet  day  pattern  was  calculated  with  the  CARP  (Calculate 
Average  Representative  Profile)  program  by  averaging  over  height  profiles  for  four  quiet  days, 
namely,  August  17-20,  2004.  The  effect  of  the  ionospheric  storm,  i.e.,  the  difference  between  the 
disturbed  and  quiet  day  data  is  quite  significant.  On  the  disturbed  day  a  series  of  electron  density 
enhancements  were  observed  between  08  UT  and  1 5  UT  at  all  stations.  These  enhancements  are 
present  over  a  large  height  interval,  from  200  km  to  450  km.  It  is  also  worth  noting  that  the 
stations  at  lower  latitudes  observe  stronger  enhancements.  This  is  in  agreement  with  the 
suggestion  by  Jones  and  Rishbeth  [1971]  that  the  intensity  of  the  positive  ionosphereic  storm 
effect  is  proportional  to  cot//  NmF2 ,  where  /  is  a  magnetic  dip  angle.  Thus,  for  the  Chilton  and 
Athens  the  ratio  between  the  expected  intensities  is  ~0.7  to  0.8.  Another  difference  between  the 
quiet  and  disturbed  days  is  the  absence  of  the  post-sunset  enhancement  of  the  electron  density 
(approximately  between  1 8  UT  and  20  UT)  on  the  storm  day.  Thus,  on  this  day,  both  positive 
and  negative  effects  are  observed  sequentially  in  the  European  sector.  It  is  also  noticeable  that 
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the  height  of  the  maximum  electron  density  moves  upward  on  the  disturbed  day  with  respect  to 
the  quiet  day  values.  To  investigate  the  height  changes  more  quantitatively  we  make  use  of  the 
foF2  and  hmF2  characteristics  deduced  from  the  electron  density  profiles.  Figure  28  illustrates 
the  comparison  of  the  layer  characteristics  recorded  on  24  August  2005  to  the  quiet  day  values. 
For  better  orientation,  the  DST  index  variation  is  plotted  in  the  top  panel  showing  storm 
commencement  and  the  following  main  and  recovery  phases. 


—  •  —  foF2 
- foF2quiet 


00  02  04  06  08  10  12  14  16  18  20  22  00 
Universal  Time 


Figure  28.  Variation  of  foF2  and  hmF2  layer  parameters  on  24  August  2005  compared  to 
the  quiet  day  averages  (black  curves).  Top  panel  shows  DST  index. 


The  data  presented  here  show  a  clear  deviation  of  the  foF2  and  hmF2  values  observed  on  the 
event  day  compared  with  the  quiet  day  pattern  emphasizing  both  positive  and  negative  (after 
16  UT)  effects  at  all  stations.  The  Ebro  station  plot  shows  the  range  of  foF2  and  hmF2  values  for 
the  four  quiet  days  used  in  calculating  the  average.  The  uplifting  observed  on  the  disturbed  day  is 
significantly  larger  than  the  height  variations  on  quiet  days.  There  are  two  major  periods  of 
strong  ionospheric  uplifting  in  these  data:  the  first  one  lasts  from  08  UT  to  14  UT  and  the  second 
one  lasts  from  1 7  UT  to  22  UT.  They  roughly  correspond  to  the  positive  and  negative  phases  of 
the  storm.  Before  looking  at  the  effect  of  ionospheric  uplifting  and  its  relation  to  the  electron 
density  enhancement/depletion  in  greater  details,  the  observations  made  in  the  American  sector 
are  shown  in  Figure  29  to  be  compared  to  the  European  sector.  Figure  29  is  in  the  same  format  as 
Figure  28  and  shows  the  foF2  and  hmF2  parameter  variations  on  the  disturbed  day  and  the 
average  quiet  values.  The  American  station  coordinates  can  also  be  found  in  Table  5.  It  is  worth 
noting  that  in  the  American  sector  the  storm  commencement  at  0613  UT  and  the  following  main 
phase  both  occur  during  local  nighttime. 
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Figure  29.  American  sector  data  for  24  August  2005  storm.  Plots  show  variations  of  foF2 
(left)  and  hmF2  (right)  as  compared  to  the  quiet  day  averages  (black  curves). 


In  Figure  29,  there  appear  to  be  times  when  foF2  and  hmF2  cannot  be  measured.  It  is 
known  that  ionospheric  storms  affect  higher  alttiude  regions  of  the  ionosphere  severely  [e.g., 
Buresova,  2005].  Thus,  during  a  negative  storm  the  foF2  parameter  can  experience  dramatic 
reductions,  while  foFl  is  only  slightly  changed.  This  often  leads  to  a  situation  in  which  foF2  < 
foFl,  and  for  an  ionosonde  this  means  that  the  F2  layer  is  “shadowed”  by  the  FI  layer.  Under 
such  circumstances,  measurements  of  both  foF2  and  hmF2  are  impossible  (this  is  referred  to  as  a 
“G-condition”  in  the  ionogram  processing  handbooks),  however,  the  upper  limit  of  foF2  can  be 
estimated  as  equal  to  foFl.  The  G-condition  is  often  observed  in  summer  daytime  at  mid¬ 
latitudes  when  the  FI  layer  is  well  developed.  Such  conditions  were  observed  on  24  August  2005 
in  the  American  sector  stations  after  sunset. 


The  behavior  of  the  ionosphere  in  the  American  sector  on  the  storm  day  is  distinctively 
different  from  that  in  the  European  one.  In  comparison  to  the  quiet  day  pattern,  we  see  an  overall 
decrease  in  the  electron  density,  i.e.,  a  negative  ionospheric  storm  was  observed  at  the  US 
stations.  In  Figure  29,  one  can  also  see  that  following  the  storm  commencement  the  stations 
observe  an  uplifting  of  the  ionosphere,  similar  to  the  one  detected  by  the  European  stations.  We 
now  need  to  investigate  this  uplifting  effect  in  both  sectors  in  more  detail  by  using  the  hmF2 
parameter  as  a  measure  of  the  ionospheric  height  variations.  Figure  30  presents  the  differences 
between  the  hmF2  values  measured  on  the  event  day  and  the  quiet  day  averages  for  both  the 
American  and  European  sectors.  These  differences  are  plotted  as  a  function  of  time  elapsed  from 
the  storm  commencement  at  0613  UT. 
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Figure  30.  Difference  between  hmF2  values  measured  on  the  event  day  and  the  quiet  day 
average  as  a  function  of  time  after  storm  commencement. 

The  approximate  beginning  of  the  F  layer  uplifting  is  indicated  by  arrows.  It  is  clear  that  in 
the  European  sector,  the  delay  between  the  storm  commencement  and  start  of  the  uplifting  is 
measureably  larger  than  that  in  the  American  sector.  The  minimal  delay  is  105  minutes  in  the 
European  sector  (Chilton)  and  less  than  60  minutes  in  the  American  sector  (Millstone  Hill). 
Moreover,  it  is  clear  (at  least  for  the  European  sector  data)  that  the  delay  increases  as  one  moves 
from  higher  to  lower  latitudes.  This  suggests  that  the  disturbance  associated  with  the  ionospheric 
storm  is  propagating  from  high  to  low  latitude.  Another  interesting  observation  comes  from  the 
American  sector  data.  Here  one  can  see  that  the  Millstone  Hill  and  Wallops  Island  have  a  time 
offset  of  one  hour  or  less,  while  the  Boulder  and  Eglin  observe  an  uplifting  at  least  2  hours  after 
the  storm  commencement.  Apparently,  longitudional  separation  plays  a  significant  role  in  the 
uplift  onset  time.  In  other  words,  the  initial  neutral  wind  surge  on  the  night  side  appears  to  be 
limited  on  the  west  side  to  a  boundary  near  01  LT  (position  of  the  Millstone  Hill  station  at  0613 
UT).  Note  that  for  this  negative  storm  in  the  American  sector,  the  uplifting  by  itself  does  not 
correlate  with  the  density  depletion  as  it  takes  another  2  to  3  hours  for  the  US  stations  to  observe 
the  foF2  decrease.  This  total  time  delay  obviously  consists  of  the  composition  bulge’s 
equatorward  progression  and  also  of  the  ionospheric  plasma  response  to  the  increased 
recombination  rate.  In  the  next  section  we  try  to  put  these  results  of  these  observations  together 
as  a  coherent  model. 
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7.3  Discussion 


Our  observations  of  the  24  August  2005  storm  are  consistent  with  the  idea  that  the  major 
driver  of  storm  effects  at  mid-latitudes  is  the  neutral  wind.  A  simplified  picture  of  the 
ionospheric  response  to  the  storm  is  as  follows  [Prolss,  1993;  Fuller-Rowell  et  al.,  1994].  A 
strong  energy  input  in  the  auroral  zone  during  the  initial  phase  of  the  storm  results  in  an  increase 
of  the  Joule  heating  that  drives  a  global  neutral  wind  surge.  Divergence  of  the  horizontal  winds 
leads  to  vertical  winds  causing  a  composition  change  in  the  upper  atmosphere,  namely,  a  change 
in  the  0/N2  ratio.  This  increase  in  O,  in  turn,  leads  to  a  higher  recombination  rate,  which  reduces 
overall  electron  density  in  the  ionosphere.  Thus,  in  the  region  of  the  composition  change,  a 
negative  ionospheric  storm  is  expected.  Such  a  region  of  increased  0/N2  ratio  is  often  referred  to 
as  a  “composition  bulge”.  This  composition  bulge  is  usually  created  on  the  nightside  in  the 
region  of  0  to  4  LT  and  then  rotates  with  the  Earth  and  is  also  pushed  by  background  neutral 
winds.  The  dayside  ionosphere  typically  observes  positive  storm  effects  as  the  horizontal 
meridional  wind  surges  originating  in  the  auroral  oval  push  the  ionosphere  upward  along  the 
magnetic  field  lines  [Rishbeth,  1998]  to  the  region  of  lower  recombination  rate  thus  leading  to 
electron  density  enhancement  after  the  start  of  the  storm.  It  is  important  to  note  that  during  the 
positive  phase  the  increase  in  the  ionozation  density  lags  the  increase  of  layer  height  by  1  to  2 
hours.  This  is  how  long  it  takes  to  build  up  an  increased  ionospheric  plasma  density  in  response 
to  the  uplifting  of  the  F  layer.  The  ionospheric  response  observed  in  our  analysis  of  the  digisonde 
data  is  in  good  agreement  with  this  explanation  of  storm  dynamics.  We  have  illustrated  this  with 
the  24  August  2005  storm  for  which  these  data  were  presented  in  detail.  First,  let  us  take  a  look 
at  the  energy  input  in  the  auroral  region.  Figure  3 1  presents  a  composite  observation  plot  made  at 
0652  UT  (storm  commencement  was  0613  UT),  showing  net  energy  flux.  This  plot  is  made 
using  the  OVATION  database  (http://sd-www.jhuapl.edu/Aurora/ovation/index.html).  We  see 
that  the  oval  was  significantly  expanded,  particularly  on  the  nightside  its  equatorward  edge 
practically  reached  the  Millstone  Hill  station.  In  the  American  sector,  all  the  digisonde  stations 
observed  a  negative  storm  effect  relatively  quickly  after  the  storm  commencement  in  agreement 
with  the  oval  observations.  Given  the  distance  from  the  northernmost  station  to  the  oval,  it  is 
natural  that  there  is  a  small  time  delay  between  the  storm  onset  and  the  observed  response  at  the 
Millstone  Hill.  This  negative  effect  is  obviously  caused  by  the  composition  bulge  surging  in 
from  the  auroral  oval.  At  the  Millstone  Hill  station,  we  also  observed  a  strong  uplifting  of  the 
ionospheric  electron  density  maximum,  which  is  evidence  of  the  presence  of  horizontal  neutral 
wind  (that  moves  the  ionosphere  upward  along  the  magnetic  field  line).  At  other  American 
stations,  the  uplifting  effect  is  not  so  clearly  detected  because  of  the  presence  of  the  G-condition 
that  made  it  impossible  to  measure  the  height  of  the  layer  peak. 
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Figure  31.  The  24  August  2005  storm.  Color  scale  shows  the  energy  input  determined  from 
radar  and  DMSP  satellite  measurement.  These  data  are  taken  from  OVATION  database. 
The  locations  of  selected  digisonde  stations  as  well  as  Bjernoya  magnetic  observatory  are 
shown  in  the  plot. 


Figure  32  presents  composite  plots  showing  the  superimposed  variations  of  AfoF2  and 
AhmF2  parameters  measured  in  the  European  sector.  We  see  that  at  all  stations  the  increase  in 
the  AfoF2  parameter  occurred  1  to  2  hours  after  the  beginning  of  the  increase  in  AhmF2 
parameter  as  is  expected  for  the  positive  storm  effect  caused  by  the  neutral  wind  disturbances.  It 
is  interesting  that  this  time  offset  is  noticeably  shorter  at  the  lower  latitudes  stations  compare  to 
those  at  the  higher  latitudes  (1  hour  at  Athens  versus  1 .5  hour  at  Chilton).  This  interesting  result 
needs  additional  research  to  be  better  understood.  Our  measurements  at  spatially  separated  points 
also  make  it  possible  to  determine  the  equatorward  propagation  speed  of  these  disturbances.  If 
measured  by  the  leading  edge,  i.e.,  by  the  beginning  of  the  increase  in  the  AhmF2  parameter, 
then  the  TAD  propagation  speed  is  between  400  m/s  and  900  m/s,  which  is  within  the  range  of 
the  expected  TAD  velocities  (for  instance,  Hocke  and  Schlegel  [1996]  report  an  average  speed 
for  the  waves  originating  in  the  auroal  latitudes  of  200-1000  m/s).  Given  the  time  offset  of  1 .45 
hour  between  the  storm  commencement  time  and  the  height  increase  at  the  Chilton  station  (see 
Figure  30)  this  TAD  propagation  velocity  puts  the  source  of  the  disturbance  causing  the  positive 
effect  in  the  European  sector  about  3000  km  north  of  the  Chilton  station.  This  corresponds  to  the 
vicinity  of  the  Bjomoya  station  in  the  middle  of  the  auroral  oval  shown  in  Figure  3 1  and  is  in 
agreement  with  the  neutral  wind-based  storm  model. 
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Figure  32.  V ariations  of  the  AhmF2  (black)  and  AfoF2  (blue)  parameters  measured  in  the 
European  sector.  The  increase  in  the  electron  density  lags  height  increases  at  all  stations.  A 
vertical  line  indicates  the  storm  commencement  time  (SST).  The  sloping  lines  connect  the 
starts  of  the  increases  in  hmF2  (black  line)  and  foF2  (blue  line). 


After  18  UT,  the  European  stations  (positioned  in  the  local  time  sector  18  to  20  LT)  observe 
the  negative  phase  of  the  storm.  Fuller-Rowell  et  al.  [1994]  have  shown  in  their  simulation  that 
1 2  hours  into  the  storm  time,  composition  disturbance  can  occupy  a  very  large  area,  essentially 
covering  the  whole  nightside  ionsophere  at  mid-latitudes.  Therefore,  it  is  not  surprising  that  so 
long  after  the  beginning  of  the  storm,  the  European  stations  have  “drifted”  into  the  zone  of  the 
disturbed  composition  and  thus  began  observing  a  negative  effect.  The  fact  that  the  easternmost 
station,  Pruhonice,  appears  to  have  detected  the  negative  phase  first  is  supporting  such  a 
scenario. 

7.4  Summary 

We  have  found  the  data  from  the  global  digisonde  network  very  useful  in  our  study  of 
ionospheric  storm  effects.  For  this  report,  we  mainly  concentrated  on  the  24  August  2005  storm. 
Digisonde  foF2  and  hmF2  data  from  the  stations  in  the  American  and  European  regions  provided 
experimental  evidence  supporting  the  storm  model  based  on  the  leading  role  of  the  neutral  wind 
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and  the  thermospheric  circulation.  The  major  results  testifying  to  the  validity  of  this  mechanism 
are  the  observation  of  the  ionospheric  uplifting  preceding  both  negative  and  positive  electron 
density  storm  phases.  Unlike  other  experimental  works  done  previously,  in  these  sounder 
observations  we  directly  measure  the  height  of  the  ionosphere  by  inverting  the  vertical  sounding 
data.  Simultaneous  observations  at  a  number  of  the  stations  is  also  a  key  part  of  our  experiment. 
The  latitudinal  chain  of  ionosondes  in  the  European  sector  observed  the  propagation  of  the 
disturbance  and  estimated  its  speed  which  places  it  into  a  class  of  the  gravity  waves/traveling 
atmospheric  disturbances.  Data  presented  here  suggest  that  the  composition  disturbance  zone 
responsible  for  the  negative  storms  has  a  tendency  to  form  around  the  sector  of  03  to  05  LT.  On 
the  other  hand,  positive  storms  prefer  to  begin  at  around  local  noontime.  Recently,  there  also  has 
been  a  lot  of  attention  paid  to  the  effects  of  the  prompt  penetration  of  the  electric  field  (e.g., 
Foster  and  Rich,  1998;  Huang  et  al.,  2005).  It  has  been  shown  that  penetrating  electric  fields  are 
also  capable  of  producing  ionospheric  storm  effects.  A  characteristic  feature  of  this  mechanism  is 
a  near-simultaneous  effect  at  a  large  range  of  latitudes.  The  effects  of  the  penetration  electric 
field  are  also  typically  short  lived,  usually  less  than  half  an  hour  (some  longer  lasting  effects 
have  been  reported  lately  as  well  [Huang  et  al.,  2005]).  The  observations  we  presented  cannot 
rule  out  the  penetration  field  mechanism,  but  the  data  collected  does  strongly  support  the  neutral 
wind  scenario  for  the  ionospheric  storm. 

8.  ROUTINE  HF  ABSORPTION  MEASUREMENBTS  USING  DIGISONDES 
8.1  Approach 

To  determine  the  diurnal  variation  of  radiowave  absorption  using  the  echo  amplitudes 
measured  with  an  ionosonde,  it  is  necessary  to  determine  the  radiated  power  of  the  specific 
sounding  system  used  for  these  measurements.  The  effective  system  gain  is  defined  as  the 
product  PtGtGr  ,  where  PT  is  the  transmitted  power  in  Watts,  Gr  and  GR  are  the  transmitting  and 

receiving  antenna  gains,  respectively.  The  received  power  is,  in  terms  of  these  factors,  given  by 
the  Friis  formula: 

p  t  r \  _  (f)  1 

(4  vfh2  L(f) 

where  k  is  the  free  space  wavelength,  h  is  the  path  length  (2  x  the  reflection  height),  and  L  is  the 
combined  losses  which  include  D-region  absorption  in  the  daytime  and  F-layer  reflection  losses. 
As  discussed  in  detail  in  the  next  section,  the  calibration  of  the  ionosonde  system  uses  the 
midnight  time  period  data  when  there  is  no  D-layer  present.  Assuming,  temporarily,  that  the 
combined  losses  L=0  during  the  selected  nighttime  hours,  the  system’s  effective  gain  at  a  given 
frequency  is: 


PtGtG 


(4;r)2  h 2 

J2 


Pr 


where  PR  is  the  received  power  derived  from  the  measured  receiver  output  voltage.  Soundings 
are  made  every  1 5  minutes  during  this  period  with  a  frequency  step  of  50  kHz.  The  frequency 
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range  for  the  calibration  presented  here  is  limited  by  the  F-layer  critical  frequency  during  the 
nighttime  period.  Here  the  frequencies  selected  for  this  analysis  are  2.5  through  5.3  MHz  in 
steps  of  0.4  MHz.  In  order  to  increase  the  number  of  data  points  at  each  frequency,  nearby 
frequencies  are  combined,  e.g.,  at  2.5  MHz,  five  frequencies  2.40,  2.45,  2.50,  2.55  and  2.60  MHz 
were  used  to  obtain  an  averaged  output  amplitude.  For  calibrating  the  sounder  system  over  the 
largest  possible  frequency  range,  the  data  for  the  months  of  May  though  August  2005  were  used. 
These  “summer”  months  had  the  highest  F-layer  critical  frequencies  during  the  nighttime 
calibration  period. 

8.2  System  Calibration 

To  calibrate  the  sounder  ionogram  amplitude  data,  the  hours  from  05  through  09  UT  (22  -  02 
LT)  were  combined  (as  well  as  the  closely  spaced  frequencies  discussed  above)  when  it  was 
reasonable  to  assume  that  the  nighttime  ionosphere  was  relatively  stable  during  these  four  hours 
on  each  of  the  nights.  In  addition,  it  was  assumed  that  there  was  no  D-region  absorption  during 
these  late  hours.  The  Digisonde  256  located  in  Boulder,  CO,  was  selected  for  this  analysis.  This 
sounder  was  chosen  for  this  preliminary  investigation  because  it  was  well  maintained  during  this 
analysis  period.  For  the  four  months  discussed  above,  the  number  of  potential  nighttime 
amplitude  measurements  exceeded  8000  per  frequency.  In  fact,  for  a  variety  of  reasons,  the 
actual  number  of  samples  used  to  determine  the  average  PTGTGR  varied  from  6500  for  2.9  MHz 
down  to  2100  for  4.1  MHz  (see  Figure  33). 


Figure  33.  Number  of  amplitude  measurements  used  to  determine  the  average  PTGTGR  as  a 
function  of  sounding  frequency. 
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For  frequencies  above  4. 1  MHz,  the  number  of  samples  over  the  four  months  became 
progressively  smaller  and  it  was  decided  at  these  higher  frequencies  to  average  nine  adjacent 
frequencies  rather  than  the  five  used  for  the  lower  frequencies.  These  small  sample  sizes  were 
caused  by  the  relatively  infrequent  occurrence  of  the  higher  F-layer  critical  frequencies.  The 
short  curve  in  Figure  33  shows  the  number  of  samples  used  at  4.5  MHz  and  above. 

Figure  34  shows  the  spread  in  the  nighttime  values  of  PTGTGR  for  2.9  MHz.  These  data  show 
an  extreme  range  of  greater  than  40  dB.  However,  93  percent  of  the  measurements  lie  between 
1 80  and  196  dB,  while  62  percent  of  the  points  lie  between  186  and  192  dB;  a  range  of  only  6 
dB.  Similar  distributions  were  calculated  for  each  of  the  other  seven  frequencies  used  in  this 
analysis. 

The  factors  that  contribute  to  this  spreading  are,  first,  focusing  and  defocusing  of  the 
reflected  signal  caused  by  irregular  structure  in  the  ionosphere,  and  second,  the  variation  in  the 
reflection  coefficient  of  the  F-layer  in  different  nights.  These  two  issues  are  addressed  in  the 
next  section. 


PTGft(cB) 


Figure  34.  Histogram  of  the  digisonde  derived  PTGTGR  values  over  the  nighttime  period 
from  May  through  August  1005  for  2.9  MHz. 
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A  similar  analysis  was  used  to  determine  PTGTGR  for  each  frequency. 

9.  FOCUSING  AND  REFLECTION  CALCULATIONS 

Addressing  the  issues  in  the  use  of  the  Friis  equation  above,  it  was  recognized  that  there  were 
two  factors  not  accounted  for  in  the  nighttime  calibration  calculation.  The  first  of  these  was  the 
focusing/defocusing  of  the  reflected  energy  by  ionospheric  irregularities  and  the  second,  the 
losses  at  the  reflection  from  the  nighttime  F-layer.  These  are  addressed  in  the  following  sections. 

9.1  Focusing  Gain 

To  better  understand  the  effects  of  ionospheric  irregularities  on  the  focusing  and  defocusing 
of  the  reflected  radio  energy,  a  program  was  written  that  generates  a  random  reflecting  surface 
with  a  specific  average  horizontal  structure  length  and  an  RMS  amplitude.  Figure  35  shows  a 
schematic  of  the  reflection  process  from  this  irregular  surface  and  illustrates  focusing  gain. 
Focusing  gain  was  defined  here  as  the  ratio  of  the  number  of  rays  reflected  from  the  irregular 

surface  that  fall  within  a  circular  area  within  a  Fresnel  radius  (  rF  =  ^JlhA  where  h  is  the 

reflection  height  and  X  is  the  free  space  wavelength)  to  the  number  of  rays  in  the  same  area 
when  the  reflecting  layer  was  smooth.  For  the  sounding  frequencies  used  in  this  study,  the 
Fresnel  radius  was  of  the  order  of  10  km.  For  these  simulations,  sounding  frequency  is  not  an 
important  factor  as  it  enters  only  in  the  Fresnel  radius  as  a  square  root  and  the  focusing  gain  ratio 
tends  to  diminish  the  effect  of  that  factor  as  well.  The  focusing  simulation  program  was 
automatically  repeated  100  times  to  develop  a  reliable  statistical  description  of  the 
focusing/defocusing  gain. 

Figure  36  shows  the  results  of  the  calculation  of  the  focusing  gain’s  100  independent  trial 
runs  (negative  gain  is  defocusing).  For  these  calculations,  the  gain  was  binned  in  3  dB  steps. 

For  the  most  shallow  amplitude  run  (RMS  amplitude  =  0.1  km),  the  result  was  close  to  what 
would  have  been  expected  from  a  smooth  reflection  surface,  i.e.,  all  the  samples  fall  in  the  0-dB 
gain  bin.  As  the  reflecting  screen  depth  was  increased  to  a  maximum  of  5  km,  the  spread  in 
focusing  gain  tends  to  increase  to  a  maximum  of  6  dB. 

The  same  program  was  then  run  (100  realizations  each)  for  four  horizontal  scale  sizes:  4.3, 
8.4,  17.0,  and  33.0  km  and  each  of  these  with  four  amplitudes:  5,  10,  15,  and  20  km.  These 
calculations  are  summarized  in  Figures  37a-37d  showing  the  percentage  of  100  realizations  with 
positive  focusing  gain  (>1.5  dB),  zero  gain  (between  -1 .5  and  +1 .5  dB)  and  negative  gain  (<-l  .5 
dB). 

The  trend  associated  with  these  calculations  is  clear.  For  each  vertical  amplitude  selected 
here,  there  is  a  horizontal  structure  size  that  yields  the  greatest  focusing  gain.  The  maximum 
gain  in  these  calculations  moves  to  larger  horizontal  structure  sizes  as  the  RMS  amplitude 
increases,  always  reaching  80%  or  more  positive  focusing  gain.  The  median  focusing  gain 
maximum  for  each  vertical  amplitude  case  is  shown  in  the  Table  6. 


48 


kkimm; 


ncrocrsiNC 


Figure  35.  Illustration  of  focusing  and  defocusing  gain  associated  with  an 
irregular  reflecting  surface. 
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Figure  36  Distribution  of  100  realizations  of  an  irregular  reflecting  surface  with  RMS 
amplitudes  varying  from  0.1  through  5.0  km.  In  all  cases  the  horizontal  scale  size  was  »  6.8 
km. 
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Figure  37.  Percentage  positive  and  negative  focusing  gain  as  a  function  of  horizontal 
structure  size  for  several  vertical  RMS  amplitudes:  a)  5  km,  b)  10  km,  c)  15  km  and  d) 
20  km. 
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Table  6.  Maximum  Focusing  Gain  and  Associated  Irregularity  Horizontal 
Structure  Size  as  a  Function  of  RMS  Amplitude 


RMS  amplitude 
(km) 

Hor.  structure  size 
for  maximum  gain 
(km) 

Max.  focusing 
gain  (dB) 

5 

8.4 

3.1 

10 

12.0 

2.8 

15 

17.0 

2.8 

20 

33.0 

3.1 

These  maximum  focusing  gains  appear  to  hold  over  a  wide  range  of  structure  sizes  for  all  the 
selected  amplitudes  and  a  good  estimate  of  overall  focusing  gain  was  3  dB. 


9.2  Nighttime  F-Layer  Reflectivity 

The  second  issue  addressed  here  is  the  assumption  made  earlier  that  the  nighttime  losses 
were  zero.  The  issue  of  reflection  from  a  layer  with  a  gradient  in  the  electron  density  is 
examined.  The  question  to  be  answered  is  how  does  the  reflection  coefficient  depend  on  the 
local  electron  density  gradient  at  the  reflection  altitude  when  there  are  no  absorbing  layers  below 
the  reflection  level? 

This  issue  was  raised  when  the  variation  of  the  nighttime  losses  was  examined  on  five 
sequential  nights  (3  May-7  May  2005).  Figure  38  shows  the  variation  of  the  losses  on  each 
individual  night  and  they  appear  to  be  relatively  constant.  However,  the  average  nighttime  loss 
varied  by  more  than  10  dB  over  these  five  nights.  This  is  well  illustrated  in  Figure  39  which 
shows  the  mean  loss  on  each  of  the  five  nights  with  bars  that  indicate  the  spread  in  the  individual 
nighttime  values. 
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Figure  38.  Boulder  measured  nighttime  loss  (05-09  UT)  on  five  sequential  days. 
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Figure  39.  Mean  loss  over  the  hours  from  05-09  UT  on  sequential  nights  from  3  May  —  7 
May  2005  (dots).  Calculated  mean  F-layer  slope  (km/MHz)  using  the  ionograms  for  the 
period  from  1  May-7  May  2005  (squares). 

Superimposed  on  the  data  curve  is  the  mean  F-layer  slope  defined  here  as  the  ratio  of  the 
height  difference  between  the  layer  peak  and  bottom  divided  by  the  measured  foF.  Although  the 
scales  for  these  two  curves  are  very  different,  the  similarity  of  the  general  shape  of  the  curves 
indicates  a  good  correlation  between  reflection  loss  and  layer  slope.  To  better  support  this 
conjecture,  a  review  of  the  theory  of  the  amplitude  of  the  reflected  waves  is  examined  in  the  next 
section. 

9.3  Reflectivity  Theory 

At  this  time,  only  a  preliminary  examination  of  this  issue  is  possible.  The  early  works  of 
Rawer  [1939],  Rawer  and  Suchy  [1967],  and  Budden  [1988]  have  investigated  the  layer 
reflection  coefficient  by  finding  solutions  to  the  wave  equation  in  a  plasma  for  the  up  going  and 
down  going  waves.  Using  a  symmetrical  Epstein  electron  density  profile,  they  calculated  the 
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reflection  coefficient  as  a  function  of  sounding  frequency  (normalized  by  the  peak  plasma 
frequency  of  the  layer)  for  several  layer  thicknesses  (see  Figure  40). 


v  -  0.02  a>p  (max)  ( ft  —  curve ) 
v  —  0. 10  cop  (max)  ( y  -  curve ) 


Figure  40.  Reflection  index  (201og|/?| )  as  a  function  of  sounding  frequency  for  normalized 

layer  thicknesses  varying  from  2  (thinnest)  to  250  (thickest)  in  three  factors  of  5.  (Rawer 
and  Suchy,  p.  194,  Figure  64). 


Here,  collision  frequencies  were  expressed  in  terms  of  its  ratio  to  the  maximum  layer  plasma 
frequency.  The  two  curves  in  Figure  40  (marked  /?.  and  y  )  were  compared  to  the  collision  free 
case  (a -curve).  Figure  41  shows  this  reflection  coefficient  as  scaled  from  the  above  curves  for 

f 

j' _  J  p, max 
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Figure  41.  Reflection  loss  (dB)  as  a  function  of  normalized  layer  thickness  for  two 
collision  frequencies. 


These  results  are  considered  preliminary  (again  the  reflection  losses  increase  rapidly  with 
increasing  layer  thickness  which  implies  a  smaller  electron  density  gradient)  so  more  effort  is 
required  to  definitively  characterize  the  nature  of  the  reflection  process  from  the  F-layer  with  a 
realistic  gradient. 
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10.  CALIBRATION  /  VALIDATION  OF  DMSP  UV  MEASUREMENTS  OF  THE  F2 
PEAK  CHARACTERISTICS  WITH  DIGISONDE  MEASUREMENTS 


The  DMSP  (Defense  Meteorological  Satellite  Program)  FI 7  weather  satellite  was  launched 
on  November  4, 2006.  It  has  two  optical  instruments  onboard: 

SSULI  -  limb  scanning  ultraviolet  imager  /  spectrometer  (built  by  the  Naval  Research 
Laboratory) 

SSUSI  -  nadir  scanning  ultraviolet  imager  /spectrometer  and  photometer  (built  by  the 
Applied  Physics  Laboratory,  Johns  Hopkins  University).  The  satellite  has  an  almost  polar  orbit 
with  a  period  of  101.9  min,  an  apogee  of  856  km,  a  perigee  of  841  km  and  an  inclination  of 
98.8°. 

The  digisonde  ground-based  electron  density  profiles  served  as  “ground  truth”  for  the 
satellite-borne  UV  optical  measurements  of  F  layer  electron  densities.  Aerospace  Corporation 
was  responsible  for  identification  of  “fly-over”  events  such  that  the  satellite  and  digisonde 
measure  the  same  region  of  the  ionosphere  at  about  the  same  time.  The  satellite  data  includes 
both  surface  and  limb-scan  measurements.  The  fly-over  calculation  results  in  at  least  one  event 
per  24-hour  period  for  each  participating  digisonde.  Each  event  represents  a  one-hour  interval  for 
which  UMLCAR  provided  manually  scaled  ionogram  data  in  SAO  4.3  format  with  15  minutes 
cadence  (can  be  different  for  non-DISS  network  sounders).  In  total,  36  stations  were  considered 
as  candidates  for  Cal/Val  project,  however,  stations  located  in  the  vicinity  of  the  South  Atlantic 
Anomaly  (SAA)  were  excluded  (6  stations).  The  following  is  the  list  of  approved  stations. 


Station  Name 

Geo. 

latitude 

Geo.  Mag. 

longitude  latitude 

Jicamarca 

-12 

283.2 

-1.78 

Kwajelain 

9 

167 

3.6 

Fortaleza 

00 

cn 

i 

322 

4.73 

Sao  Luis 

-2.6 

315 

6.41 

Hainan  Island 

19.4 

109 

9.14 

Chung-Li 

25 

121.2 

14.99 

Madimbo 

-22.38 

30.88 

-19.79 

Wuhan 

30.6 

114.4 

20.4 

Osan 

37.1 

127 

27.33 

Ramey 

18.5 

292.9 

28.72 

Boulder 

40 

254.7 

-31.25 

Learmonth 

-21.3 

114 

-31.99 

Athens 

38 

23.6 

36.35 
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San  Vito 

40.6 

17.8 

39.89 

Eglin 

30.4 

273.2 

40  2(T 

El  Arerosillo 

37.1 

353.333 

40.8 

Point  Arguello 

34.7 

239.4 

41.04 

Dyess 

32.5 

260.3 

41.41 

Rome 

41.9 

12.5 

42.09 

Irkutsk 

52.4 

104 

42.17 

Roquetes/T  ortosa 

40.8 

0.3 

43.19 

Bundoora 

-37.72 

145.05 

-45.62 

Pruhonice 

50 

14.6 

49.55 

Dourbes 

50.1 

4.6 

51.43 

Millstone  Hill 

42.6 

288.5 

52.86 

Chilton 

51.5 

359 

53.76 

Wallops 

37.9 

284.5 

54.04 

Fairford 

51.7 

358.5 

54.04 

Juliusruh 

54.6 

13.4 

54.17 

Goose  Bay 

53.3 

300 

63.27 

Requests  for  following  stations  have  been  processed. 


Station  Name 

Geo. 

latitude 

Geo.  Mag. 

longitude  latitude 

Jicamarca 

-12 

283.2 

-1.78 

Kwajelain 

9 

167 

3.6 

Fortaleza 

-3.8 

322 

4.73 

Sao  Luis 

-2.6 

315 

6.41 

Madimbo 

-22.38 

30.88 

-19.79 

Ramey 

18.5 

292.9 

28.72 
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Table  7  shows  the  status  of  the  data  requests  for  the  period  of  January  1,  2007,  until  September 
30,  2007. 


Table  7.  CalVal  Digisonde  Data  Request  Status 


Request  status 

Number  since  January  1,  2007 

Total  requests 

138 

Requests  loaded  to  DIDBase 

122 

Requests  not  loaded  (ADRES  system  is 
waiting  for  data  availability) 

16 

Not  available  (no  data)  requests 

14 

Loaded  but  not  scaled  requests 

0 

Scaled  and  reported  requests 

122 

Scaled  and  reported  ionograms 

533 
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